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Abstract  of  Dissertation  Presented  to  the 
Graduate  Council  of  the  University  of  Florida  in  Partial 
Fulfillment  of  the  Requirements  for  the  Degree  of  Doctor  of  Philosophy 

PLASTIC  STRESS  WAVES  IN  PRESTRESSED  THIN-WALLED  TUBES: 

NUMERICAL  ANALYSIS  BY  RATE-DEPENDENT  THEORIES 

By 

Arun  K.  Banerjee 
December,  1972 

Chairman:  Dr.  L.  E.  Malvern 

Major  Department:  Engineering  Science,  Mechanics, 

and  Aerospace  Engineering 

Numerical  solutions  are  reported  for  plastic  stress  wave 
propagation  in  prestressed  thin-walled  tubes  due  to  torsional  impact 
and  combined  longitudinal  and  torsional  impact.  Some  specific  unload- 
ing problems  associated  with  combined  stress  waves  are  also  considered. 
All  the  solutions  are  based  on  rate-dependent  constitutive  equations. 
Radial  inertia  has  been  neglected  throughout. 

The  torsional  loading  wave  solutions  are  compared  with  the 
experimental  data  of  Yew  and  Richardson  (1969).  Two  semilinear  consti- 
tutive equations  incorporating  linear  and  exponential  overstress  and 
one  quasilinear  constitutive  equation  are  examined.  An  integro- 
differential  approach  is  used  for  the  exponential  overstress  problem, 
and  a higher  order  Courant,  Isaacson,  Rees  method  for  the  quasilinear 
model.  Incremental  wave  speed  and  large  level  strain  data  are  better 
matched  by  all  the  three  rate-dependent  solutions — the  quasilinear 
model  giving  the  best  overall  agreement — than  a rate-independent  solu- 
tion based  on  the  quasistatic  stress-strain  relation. 
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The  combined -stress  loading  wave  -solutions  are  compared  with 
the  experimental  results  and  rate-independent  solution  given  by  Lipkin 
and  Clifton  (1970).  The  plastic  strain-rate  equations  are  derived 
from  a viscoplastic  potential  function  as  suggested  by  Rice  (1970). 

A generalized  combined  hardening  formulation,  based  on  the  work  of 
Goel  and  Malvern  (1970) , is  given.  Numerical  solutions  are  obtained 
by  a Newton-Raphson  scheme  employing  Gauss-Jordan  elimination.  Theo- 
retical solutions  for  the  special  case  of  isotropic  hardening  retain 
the  good  features  of  the  rate-independent  solution  while  bringing  out 
qualitatively  the  non-constant  state  response  in  the  experiments  where 
the  rate-independent  theory  predicted  constant  states.  While  the  prop- 
agation speed  for  longitudinal  strains  predicted  by  both  the  theories 
disagrees  with  the  experimental  data,  the  propagation  speed  of  large 
level  incremental  shear  strains  reported  from  the  experiments  is  better 
matched  by  the  rate-dependent  solution  given.  Some  preliminary  results 
based  on  kinematic  hardening  and  combined  isotropic  and  kinematic 
hardening  are  also  presented. 

Some  specific  unloading  problems  concerning  combined  stress 
waves  are  studied  next.  It  is  shown  that  the  initial  slope  and  the 
instant  of  initiation  of  the  loading/unloading  boundary  based  on  a 
rate-dependent  theory  is  generally  different  from  that  given  by  a 
rate-independent  theory,  unless  the  relaxation  boundary  has  already 
been  reached  when  decrease  of  stress  occurs  at  the  boundary.  Some 
illustrative  mixed  boundary  value  problems^ are  solved  showing  unloading 
in  a prestressed  semiinfinite  tube  and  a tube  of  finite  length. 

Numerical  solutions  bring  out  some  interesting  details  that  are  not 
easily  visualized  in  advance. 


xiii 


CHAPTER  1 


INTRODUCTION 


1.1  General  Background 

Dynamic  plasticity  is  the  study  of  the  effect  on  Solid  bodies 
of  momentary  loads  beyond  the  elastic  limit.  Examples  of  such  loads 
on  structures  are  those  caused  by  earthquakes,  explosive  blasts,  or 
collisions  between  moving  bodies.  The  subject  is  relatively  new,  with 
a history  of  about  30  years,  but  within  this  period  it  has  experienced 
a most  vigorous  growth.  It  has  been  the  topic  for  one  international 
symposium  (Kolsky  and  Prager,  1963),  numerous  conferences  in  several 
countries,  one  textbook  (Cristescu,  1967)  with  about  800  references, 
and  various  survey  papers.  In  the  following,  an  attempt  is  made  to 
recount  the  most  prominent  phases  of  development  of  the  field,  while 
for  fuller  surveys  the  reader  should  see  the  review  articles  of 
Hopkins  (1961) , Craggs  (1961) , Goldsmith  (1961) , and  Cristescu  (1968) . 

The  first  systematic  attempt  at  a serious  study  of  the  problem 
was  made  by  Donnell  (1930).  He  considered  the  propagation  of  a uni- 
axial stress  pulse  in  a thin  bar  exhibiting  a bilinear  stress-strain 
relationship.  One  significant  finding  was  that  the  nonlinearity  was 
responsible  for  a change  in  the  shape  of  the  stress  pulse  along  the 
length  of  the  bar.  Further  investigations  were  not  to  take  place 
until  a decade  later. 
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World  War  II  emphasized  the  need  for  vigorous  investigation  in 
the  general  area  of  impact  and  the  consequent  damage  on  structures. 

The  basic  problem  of  longitudinal  impact  of  thin  rods  was  studied  by 
Taylor  (1940)  in  the  United  Kingdom,  von  Karman  (1942)  in  the  United 
States,  and  Rakhmatulin  (1945)  in  the  Union  of  Soviet  Socialist  Repub- 
lics, working  independently  of  each  other.  This  has  been  recognized 
as  marking  the  true  beginning  of  the  study  of  the  propagation  of  plastic 
stress  waves  in  solids.  A finite  stress-strain  relation  was  used  in 
these  analyses  with  the  assumption  that  it  did  not  vary  with  the  rate 
of  straining.  The  main  results  were  as  follows.  A disturbance  at 
the  boundary,  such  as  a constant  velocity  impact,  propagates  itself  by 
stress  waves . the  leading  front  of  which  moves  with  elastic  speed,  and 
immediately  at  this  wave  front  there  is  a discontinuity  in  strain  corre- 
sponding to  the  initial  yield  strain.  Behind  this  wave  there  is  a 
continuous  wave  in  which  the  velocity  of  propagation  diminishes  from 
elastic  to  plastic  wave  speed  corresponding  to  the  maximum  strain 
developed,  as  the  strain  increases  from  the  initial  yield  value  to  the 
maximum  strain  reached  (which  is  a function  of  impact  velocity  and  the 
stress-strain  relation).  The  strain  then  remains  constant  forming  a 
"plateau"  up  to  the  impact  end  of  the  bar. 

Experimental  verification  of  the  above  theory  was  carried  out 
by  Duwez  and  co-workers  (1942,1947).  In  general,  the  prediction  of 
a plateau  of  strain  and  the  relation  between  the  impact  velocity  and 
the  maximum  plastic  strain  reached  were  fairly  well  verified  by  the 
experiments.  The  only  significant  and  systematic  discrepancy  observed 
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was  in  the  distribution  of  strain.  In  the  tension  impact  tests  the 
maximum  residual  strain  was  smaller  than  predicted  by  the  theory,  and 
the  observed  stress-time  variation  at  the  fixed  end  during  impact 
showed  the  stress  there  was  greater  than  the  theory  predicted.  It  was 
suggested  that  these  discrepancies  were  due  to  the  use  in  the  theory 
of  the  invariant  relation  between  stress  and  strain,  independent  of 
strain  rate.  Historically  speaking,  this  is  the  background  for  the 
subsequent  development  of  theories  that  included  the  strain-rate  effect 
as  well. 

1. 2 Rate  Dependence:  Theory  and  Experiment 

The  important  first  step  in  the  development  of  rate-dependent 
theories  was  made  by  Sokolovsky  (1948)  and  by  Malvern  (1949,1951), 
independently  of  each  other.  Sokolovsky's  theory  assumed  perfect 
plasticity,  while  Malvern  considered  the  more  general,  work  hardening 
solid.  The  constitutive  law  proposed  by  Malvern  was  of  the  elastic- 
visco-plastic  form  given  by 

E*e  = b + g(a,  e) 

where  g is  an  arbitrary  function  expressing  the  strain-rate  dependency. 
A special  case  was  solved,  choosing  g(cr,e)  = k(a-f(e)),  where  a-f(e) 
was  called  the  "overstress,"  that  is,  the  excess  of  dynamic  stress 
over  the  static  yield  stress.  In  these  relations,  elastic  strain  is 
assumed  independent  of  the  strain  rate  an£  no  account  is  taken  of 
lateral  inertia.  Malvern's  theory  has  certain  immediate  implications. 
First,  there  is  an  apparent  increase  in  the  value  of  the  initial  yield 


4 


stress  at  high  strain  rates  because  plastic  strain  does  not  occur 
instantaneously  but  takes  time  in  which  to  become  appreciable. 

Secondly,  small  increments  of  strain  superposed  on  a static  plastic 
prestrain  are  propagated  with  elastic  wave  velocity  and  not  with  the 
prestrain  plastic  wave  velocity  predicted  by  the  strain-rate  inde- 
pendent theory  of  Taylor-von  Karman-Rakhmatulin. 

Regarding  the  first  point,  apparent  increase  of  yield  stress, 
the  experiments  of  Campbell  (1954)  with  steel  specimens  showed  that 
the  stress  necessary  to  cause  yield  is  about  twice  as  great  under  the 
conditions  of  an  impact  test  as  that  required  under  normal  static  con- 
ditions. Repeated  dynamic  loading  tests  on  mild  steel,  conducted  by 
Taylor  (1957) , showed  that  the  dynamic  yield  stress  increases  with  the 
rate  of  strain.  More  recently,  Frantz  and  Duffy  (1971)  have  shown, 
from  tests  on  aluminum  giving  a superposed  high  strain  rate  on  a "static" 
rate  and  a static-plastic  prestrain,  that  the  initial  response  to  the 
strain-rate  increment  is  elastic,  followed  by  yielding  behavior  remi- 
niscent in  appearance  to  an  upper  yield  stress. 

The  velocity  of  propagation  of  incremental  waves  in  longitud- 
inal impact  of  a prestressed  steel  bar  was  measured  first  by  Bell  (1951) 
and  by  Bell  and  Stein  (1962) , and  was  found  to  be  equal  to  the  elastic 
velocity.  Similar  conclusions  were  arrived  at  for  incremental  waves 
in  copper  by  Sternglass  and  Stuart  (1953)  and  for  lead  by  Alter  and 
Curtis  (1956).  These  results  thus  agree  with  the  predictions  of 
Malvern's  theory.  However,  Bell  and  Stein  (1962)  attributed  the 
identical  observed  effect  in  annealed  aluminum  to  serrations 
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(Portevin  - le  Chatelier  effect)  in  stress-strain  curves  and  not  to 
strain-rate  effects.  Bell  (1968)  has  in  fact  summarized  the  results 
of  a series  of  tests  on  constant  velocity  impact  of  metal  bars  in  free 
flight.  For  certain  metals,  including  aluminum  in  the  fully  annealed 
condition,  the  velocities  of  propagation  of  given  levels  of  strain  were 
found  to  be  independent  of  strain  rate,  and  the  particle  velocities 
satisfied  a certain  integral  relationship  arising  in  the  rate-independ- 
ent theory  between  velocity,  wave  speed,  and  the  finite  strain.  Thus 
the  conclusion  was  that  a single  dynamic  stress-strain  relationship 
exists  so  that  a strain-rate  independent  theory  is  applicable.  This  is 
the  basis  of  the  well-known  Bell  (1963)  parabola.  However,  Bell  main- 
tains that  the  dynamic  stress-strain  relation  is  not  necessarily  the 
same  as  the  quasistatic  relation  and  thus,  as  Hopkins  (1961)  observes, 
although  there  is  strain-rate  independency  for  the  dynamic  experimental 
conditions  themselves,  there  is  strain-rate  dependency  for  these  condi- 
tions in  comparison  with  the  quasistatic  ones. 

Hauser,  Simmons  and  Dorn  (1961)  studied  the  strain-rate  effects 
in  plastic  wave  propagation  in  aluminum,  using  a split  Hopkinson  pres- 
sure bar,  and  concluded  that  the  Malvern  theory  is  a good  approximation 
for  the  material  constitutive  equation,  although  the  function  g(<j,S) 
may  not  be  a simple  function  of  the  overstress,  c-f(e).  Lindholm 
(1964)  showed  from  similar  tests  that  the  prior  strain-rate  history  of 
the  specimen  had  a significant  effect  on  plastic  flow  behavior  when 
reloading  was  at  a high  strain  rate.  Also,  dynamically  reloaded  speci- 
mens indicated  an  annealing  recovery  effect  with  a characteristic  time 
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of  the  order  of  seconds.  These  findings  do  not  agree  with  the  concept 
of  a single  dynamic  stress-strain  relation.  The  strain-rate  effect 
also  was  positively  substantiated  by  the  numerous  experiments  of 
Ripperger  and  his  associates  (Ripperger  (1960) , Tapley  and  Plass  (1961) , 
Tapley  (1962) , Karnes  (1963) , Ripperger  and  Watson  (1967)  , Baker  and 
Yew  (1966)).  The  experimental  work  of  Bodner  and  Symonds  (1962)  on 
impact  loading  of  cantilever  beams  of  mild  steel  and  aluminum  also 
showed  that  better  agreement  between  experiment  and  theory  could  be 
obtained  by  incorporating  in  the  latter  an  empirical  account  of  the 
effect  of  strain  rate  on  the  relation  between  bending  moment  and  curva- 
ture. 

The  controversy  over  rate-dependence  rages  still.  While  Bell 
(1968),  from  his  experiments,  sees  an  overwhelming  evidence  against 
rate  effect,  other  experimenters  do  feel  the  need  for  including  rate 
sensitivity  for  interpreting  their  data.  It  has  been  generally  accepted 
that  while  the  rate-independent  theory  is  adequate  for  explaining 
finite  strain  data,  the  rate-dependent  theories  have  the  capability 
of  explaining  incremental  strain  response.  At  least  on  one  particular 
aspect,  however,  the  existence  of  a strain  plateau,  the  two  viewpoints 
give  the  same  result.  While  it  was  initially  believed  that  the  Malvern 
theory  cannot  predict  a plateau,  Rubin  (1954)  showed  that  it  does 
explain,  only  by  asymptotically  approaching,  the  rate-independent 
plateau.  Later  work  by  Wood  and  Phillips  (1967)  and  Efron  (1964) 
confirmed  this.  A more  recent  theoretical*  work  by  Suliciu,  Malvern 
and  Cristescu  (1972)  has  placed  restrictions  on  the  form  of  the 
constitutive  equations  under  which  a plateau  can  occur  and  has  shown 
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that  with  the  usual  form  of  the  Malvern-type  semilinear  models  only 
an  asymptotic  plateau  is  possible. 

A theory  lying  between  the  two  extremes  of  rate-independent 
theory  and  Malvern's  theory  is  the  one  proposed  by  Cristescu  (1963). 
Lubliner  (1964)  also  used  a similar  equation.  It  has  the  form 

e = (i  + cp(CT,e))cr  + ijj(cr,e) 

£i 

with  the  implication  that  a certain  part  of  the  plastic  strain  is 
also  developed  instantaneously.  Working  on  Bell's  experimental  data, 
Cristescu  (1972)  has  shown  that  the  above  constitutive  equation,  with 
special  forms  for  the  functions  cp  and  \Jr , can  fit  data  on  various 
aspects  of  the  problem  of  longitudinal  impact  between  two  thin  rods. 

In  the  present  dissertation,  the  analyses  will  be  based  on 
the  rate-dependent  theories  of  Malvern  and  Cristescu,  suitably  modi- 
fied to  problems  in  one  or  two  stress  components. 

1. 3 Previous  Work  Related  to  the 
Problems  Considered  in  Thesis 

Previous  investigations  on  plastic  loading  waves  in  torsion 
and  combined  stress  and  the  problems  of  unloading  associated  with 
combined  stress  will  be  considered  below. 

Torsional  Waves.  Experiments  in  torsional  waves  were 
undertaken  to  eliminate  any  complicating  effect  due  to  radial 
inertia  that  arises  in  longitudinal  impact.  Baker  and  Yew  (1966) 
developed  a torsional  split  Hopkinson  bar  and  compared  their 
experimental  shear  strain-time  records  with  both  rate-independent 
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theory,  using  the  quasistatic  stress-strain  relation  and  rate-dependent 
theory  with  linear  overstress  as  used  by  Malvern.  The  conclusion  was 
that  the  rate-dependent  theory  gave  better  agreement. 

Convery  and  Pugh  (1968) , at  the  suggestion  of  Craggs  (1961) , 
studied  incremental  waves  in  torsion  by  a new  method  in  which  a tubular 
test  specimen  together  with  a concentric  bar  of  a brittle  material 
was  twisted  slowly  such  that  when  the  specimen  was  stressed  beyond  its 
yield  the  brittle  bar  broke  suddenly  and  transmitted  a plastic  tor- 
sional stress  increment  along  the  speciman.  The  speed  of  propagation 
of  this  incremental  stress,  in  both  copper  and  mild  steel,  was  found 
to  be  the  same  as  the  elastic  shear  wave  velocity.  This  is  defin- 
itely a very  convincing  proof  of  the  strain-rate  effects  in  these 
metals,  justifying  the  use  of  rate-dependent  theories. 

Yew  and  Richardson  (1969)  also  studied  incremental  plastic 
waves  in  torsion  with  copper  specimens,  but  as  their  experimental  data 
forms  the  basis  of  the  theoretical  analyses  reported  in  Chapter  2,  it 
will  be  discussed  in  that  context. 

Combined  Stress  Waves.  To  the  best  of  this  author's  knowledge, 
all  of  the  theoretical  works  published  so  far,  with  the  exceptions 
of  those  by  Perzyna  (1962)  and  Hsu  (1972)  to  be  described  later,  are 
based  on  the  inviscid  (that  is,  rate-independent)  plasticity 
theory. 

Craggs  (1957)  first  gave  an  analysis,  assuming  small  strain 
and  rotation,  of  combined  stress  plane  waves  in  an  unbounded  medium 
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and  showed  that,  quite  generally,  there  are  two  wave  speeds,  a fast 

wave  speed  c and  a slow  wave  speed  c , such  that 
f s 


c 

s 


< c„ 


< c. 


< c 
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where  c is  the  elastic  dilatational  wave  speed  and  c is  the  elastic 
1 ^ 

shear  wave  speed. 

Rakhmatulin  (1958)  studied  the  problem  of  edge  impact  of  two 
plates  for  which  the  directions  of  the  velocities  of  the  two  plates 
are  in  the  plane  of  the  plates  and  are  oblique  to  the  impacting  faces. 
His  analysis,  based  on  Hencky's  total  deformation  theory,  showed  that 
a uniaxial  stress  state  would  propagate  first  as  a simple  wave,  fol- 
lowed by  a constant  state;  then  a combined  stress  state  would  propa- 
gate as  a strong  discontinuity,  followed  by  a constant  state. 

Cristescu  (1959),  using  the  same  equations  as  Rakhmatulin  and 
assuming  a material  for  which  the  slope  of  the  stress-strain  curve  in 
shear  is  less  than  the  secant  modulus,  found  two  types  of  combined 
stress  waves.  He  showed  that  a given  level  of  compressive  strain 
propagates  faster  in  the  faster  of  these  two  combined  stress  waves 
than  it  would  in  a plastic  wave  for  a uniaxial  stress  state;  similarly, 
in  the  slower  combined  stress  wave  the  velocity  of  propagation  of  a 
given  level  of  shear  strain  was  shown  to  be  faster  than  in  a plastic 
wave  of  pure  shear. 

Bleich  and  Nelson  (1966)  gave  a closed  form  solution  to  the 
problem  of  a half  space  with  a uniformly  distributed  step  load  of 
pressure  and  shear,  for  an  elastic-perf ectly  plastic  material. 
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Clifton  (1966)  presented  the  solution  for  a thin  walled  semi- 
infinite tube  with  combined  longitudinal  and  torsional  step  loading 
for  an  elastic-plastic  material  with  isotropic  hardening.  For  various 
initial  and  boundary  conditions  regarding  the  stress  state,  he  showed 
the  existence  of  simple  wave  regions  (where  the  state  variables  prop- 
agate with  constant  velocity)  in  the  characteristic  plane  with  asso- 
ciated slow  or  fast  wave  speeds  mentioned  in  the  beginning. . In  con- 
trast with  the  results  of  Rakhmatulin  and  Cristescu,  this  theory  pre- 
dicted a wave  of  strong  discontinuity  occurring  only  in  elastic  regions, 
and  also  a uniaxial  plastic  wave  propagating  as  the  leading  wave  due 
to  a step  load  of  normal  stress  and  shear. 

Ting  and  Nan  (1969)  generalized  the  half-space  problems  of 
Bleich  and  Nelson  for  an  isotropic  hardening  material.  Ting  (May, 

1968)  also  solved  a boundary  value  problem  for  a half  space  with  two 
shear  loadings  only,  and  in  another  paper  (Dec.,  1968)  gave  a formula- 
tion of  the  partial  differential  equations  for  a general  plane  wave 
in  a half  space.  Nan  (1968)  obtained  numerical  solutions  for  a half 
space  subjected  to  two  impact  loads  in  shear  or  one  compression  and 
one  shear. 

Goel  and  Malvern  (1970)  considered  the  problem  of  Clifton,  but 
introduced  in  their  analysis  of  wave  propagation  the  concept  of  com- 
bined kinematic,  and  isotropic  hardening  for  a material  characterized 
by  a nonlinear  stress-strain  curve.  Simple  wave  stress  trajectory 
calculations  showed  that  the  hardening  law  can  make  a significant 
difference  in  the  results.  The  same  authors  (1971)  considered  in 
another  paper  the  propagation  of  elastic/plastic  plane  waves  with 
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combined  compressive  and  two  shear  stresses  in  a half  space  of  an 
isotropic  material.  Simple  wave  solutions  and  an  example  of  a non- 
centered  simple  wave  were  presented. 

Lipkin  and  Clifton  (1970)  made  an  experimental  study  and 
rate-independent  analysis  of  combined  stress  waves  in  a pretorqued 
thin  tube.  This  experiment  forms  the  basis  for  the  theoretical  work 
reported  in  Chapter  3 and  will  be  discussed  there. 

A rate-dependent  analysis  of  combined  stress  waves  was  made  by 
Perzyna.  In  a sequence  of  reports  (April,  May,  June,  1962),  he  first 
generalized  the  one-dimensional  constitutive  equation  for  rate  sensi- 
tive plastic  materials  to  general  states  of  states  and  then  showed  that 
the  problems  of  spherical  waves,  cylindrical  radial  waves,  cylindrical 
shear  waves  and  plane  waves  in  a half  space  could  be  reduced  to  the 
same  mathematical  formulation.  Methods  of  solution  were  also  discussed. 

Combined  Stress  Problems  with  Unloading.  No  published  work, 
except  that  of  Hsu  (1972) , is  available  that  considers  combined  stress 
unloading  in  a rate  sensitive  medium.  The  following  papers  are  based 
on  classical  plasticity  theory.  Clifton  (1968)  considered  the  problem 
of  elastic/plastic  boundaries  in  combined  longitudinal  and  torsional 
plastic  wave  propagation.  Equations  giving  the  restrictions  on  the 
possible  speeds  of  elastic/plastic  boundaries  were  obtained  and  these 
restrictions  were  shown  to  depend  on  the  type  of  discontinuity  at  the 
boundary  and  on  whether  loading  or  unloading  was  occurring.  The  gen- 
eral features  of  the  discontinuities  associated  with  loading  and 
unloading  boundaries  were  established  and  examples  were  presented  of 
unloading  boundaries  overtaking  simple  waves. 
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Ting  (1969)  considered  the  initial  slope  problem  of  elastic/ 
plastic  boundaries  in  combined  longitudinal  and  torsional  wave  propa- 
gation. The  stresses  applied  at  the  boundary  were  assumed  to  be  con- 
tinuous and  the  initial  slope  of  the  boundaries  was  determined  analy- 
tically for  all  possible  combinations  of  the  time  derivatives  of 
stresses  before  and  after  the  elastic/plastic  boundaries  occurred. 
Among  other  interesting  results,  it  was  shown  that  more  than  one 
elastic  and  plastic  region  could  be  generated  near  the  boundary  even 
though  the  stresses  at  the  boundary  were  continuously  increasing  to 
a higher  yield  surface. 

In  a more  recent  paper,  Ting  (1972)  considered  the  solution 
near  the  impact  end  of  a prestressed  tube  subjected  to  a discontin- 
uous combined  longitudinal  and  torsional  loading  at  time  t = 0.  The 
loads  after  t = 0 were  assumed  continuous.  Solutions  were  obtained 
for  all  possible  combinations  of  the  discontinuous  loadings  at  t = 0 
and  all  possible  variations  after  t = 0.  The  characteristics  of  the 
solutions  in  each  region,  the  initial  speed  and  in  some  cases  the 
initial  acceleration  of  the  boundary  between  two  regions  were  given. 

Hsu  (May,  1972)  conducted  tests  on  combined  stress  waves  in 
a pretorqued  tube  of  alpha-titanium, which  was  found  to  be  strongly 
rate  sensitive.  Numerical  solutions,  involving  unloading  from  impact 
end  were  carried  out,  based  on  a visco-plastic  model  related  to  an 
elementary  dislocation  mechanism,  and  substantial  agreement  was  found 


with  the  experiment. 
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1. 4 Outline  of  the  Present  Study 

The  study  reported  in  this  thesis  represents  a numerical 
investigation  of  incremental  loading  waves  in  torsion  and  combined 
stress  and  some  unloading  problems  associated  with  the  latter.  In  all 
cases  the  purpose  was  to  know  the  consequences  of  various  rate-dependent 
constitutive  laws  with  semilinear  or  quasilinear  formulation  and  spe- 
cific hardening  assumptions.  As  a guide  as  well  as  a check  to  the  theo- 
retical results,  an  effort  was  made  to  match  experimental  data  avail- 
able for  the  problems  considered,  even  though  the  materials  used  in  the 
experiments  were  reportedly  not  very  rate  sensitive. 

In  Chapter  2,  a theoretical  analysis  is  given  to  interpret  the 
experimental  data  on  incremental  waves  in  torsion  as  reported  by  Yew 
and  Richardson  (1967,1969).  Three  different  constitutive  equations  are 
examined.  The  first  one,  a semiliiiear  form  using  linear  overstress,  is 
readily  integrated..  The  second  one,  also  semilinear  but  involving  an 
exponential  function  of  the  overstress,  gives  rise  to  serious  instabil- 
ity problems  in  computation.  The  required  stability  is  built  into  the 
computations  by  first  giving  a variation  of  parameter  solutions  to 
the  constitutive  equation  and  then  introducing  it  into  the  wave  equa- 
tion to  obtain  an  integro-dif f erential  equation,  which  was  solved  by 
the  method  of  characteristics.  The  third  constitutive  law  considered 
is  of  the  Cristescu  (1972)  type.  Solutions  to  the  resulting  system  of 
quasilinear  equations  were  obtained  by  the  second-order  accurate 
Courant,  Isaacson,  Rees  method  developed  by  Ranganath  and  Clifton 


(1971). 
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In  all  three  cases,  stress  at  the  boundary  (boundary  strain  was 

the  input)  and  incremental  strains  at  two  other  stations  were  computed. 
Incremental  strain-time  records  from  the  experiments  were  compared  with 
the  theoretical  solutions.  A reasonable  agreement  was  obtained  for 
high  strain  levels,  but  not  so  for  the  small  strain  ones.  It  is  shown 
that  such  strain-time  records  or  velocity  of  propagation  measurements 
obtained  from  them,  by  themselves,  are  not  sufficient  for  determining 
the  constitutive  relation. 

In  Chapter  3,  an  attempt  is  made  to  explain  the  experimental 
results  on  combined  stress  waves  reported  by  Lipkin  and  Clifton  (1970) 
from  a rate-dependent  point  of  view.  A Perzyna  (1962) -type  constitu- 
tive equation  is  derived  from  an  appropriate  viscoplastic  potential 
function  as  postulated  by  Rice  (1970).  A formulation  is  given  for 
combined  isotropic  and  kinematic  hardening.  The  governing  system  of 
nine  quasilinear  partial  differential  equations  is  integrated  by  the 
method  of  characteristics,  utilizing  a Newton-Raphson  scheme. 

Numerical  results  are  presented  for  three  cases  of  hardening 
mechanisms — isotropic,  kinematic  and  combined.  Solutions  for  normal 
and  shear  stresses,  (a-T)  coordinates  of  the  center  of  the  yield 
surface,  the  size  parameter  describing  the  yield  surface,  the  equiv- 
alent overstress,  normal  strain,  and  the  incremental  shear  strain  are 
reported.  Experimental  records  of  compressive  strain  and  incremental 
shear  strain  reported  at  one  station  are  in  better  agreement  with 
present  theory  using  isotropic  hardening  than  the  rate  independent 
solutions  given  by  Lipkin  and  Clifton.  Convergence  difficulties  of 
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the  theoretical  solution  arising  out  of  an  attempt  to  describe  a not 
very  rate  sensitive  material  behavior  by  a rate-dependent  theory  seem 
to  militate  against  a better  agreement  with  the  experiments  in  the 
other  two  cases  of  hardening. 

In  Chapter  4,  some  unloading  problems  in  combined  stress  waves 
are  considered.  First  the  initial  slope  problem  for  different  shapes 
of  the  boundary  stress  pulses  is  studied  and  it  is  shown  that  for  a 
rate  sensitive  medium  the  loading/unloading  boundary  may  or  may  not 
start  at  the  instant  of  unloading.  It  is  shown  that  for  the  stress 
conditions  that  produce  unloading  in  a rate  insensitive  material 
(Ting,  1969),  initial  slope  behavior  of  the  loading/unloading  boundary 
in  a rate  sensitive  material  is  generally  identical  with  the  rate 
insensitive  situation  only  if  a sufficient  dwell  time  is  allowed  for 
the  associated  relaxation  phenomena  to  take  place  before  unloading 
occurs. 

Next,  a mixed  boundary  value  problem  is  solved  for  a semi- 
infinite tube  with  velocities  specified  during  loading  and  stresses 
during  unloading  at  the  boundary.  Work  hardening  behavior,  overstress 
data  at  the  boundary,  normal  strain  and  incremental  shear  strain  at 
two  stations  are  given  for  two  pulse  shapes. 

Finally,  the  same  problem  is  considered  for  a finite  tube  with 
the  other  end  constrained  against  rotation  but  free  for  axial  motion. 
Isolated  domains  of  loading  and  unloading  are  obtained.  Some  inter- 
esting and  not  intuitively  expected  result*  at  the  two  boundaries  and 
the  two  interior  stations  are  reported. 

In  Chapter  5 the  results  are  summarized  and  comments  and 


conclusions  are  given. 


CHAPTER  2 


INCREMENTAL  STRAIN  WAVES  IN  TORSIONAL  IMPACT 

The  objective  of  the  material  presented  in  this  chapter  is 
to  analyze,  from  the  point  of  view  of  strain-rate  dependent  theory, 
results  of  an  experiment  on  incremental  plastic  stress  waves  in  torsion 
as  reported  in  the  literature.  In  particular,  the  consequences  of  vari- 
ous constitutive  laws  on  the  theoretical  solutions  concerning  wave  prop- 
agation will  be  studied  and  compared  with  the  experimental  data. 

2.1  The  Experiments  of  Yew  and  Richardson 

The  problem  of  Yew  and  Richardson  (1967,1969),  who  describe  the 
details  of  their  experimental  arrangement,  is  considered.  An  annealed 
copper  tube,  0.500  in.  O.D.  and  0.433  in.  I . D.  and  24  in.  long,  is 
first  cold-worked  to  a state  of  prestress  exceeding  the  yield  stress 
by  subjecting  it  to  a very  slowly  increasing  torque  (straining 
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rate  3 x 10  /sec).  Next  an  impulsive  torque  is  generated  in  the  tube 

by  means  of  the  torsional  impact  machine  developed  by  these  authors. 

Electrical  wire  resistance  strain  gauges  mounted  on  the  specimen  tube 

at  three  locations,  0.38  in.,  1.50  in.  and  2.75  in.,  respectively,  from 

the  impact  end,  monitor  on  an  oscilloscope  screen  the  incremental  shear 
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strain  against  time.  For  an  elastic  shear  wave  speed  of  8.3  x 10 
in. /sec,  as  reported  by  these  authors  (1967),  the  first  unloading  wave 
reflected  from  the  free  end  of  the  tube  takes  545  microsec  to  reach 
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the  gauge  farthest  from  the  impact  end.  Thus,  for  the  duration  of 

impact  for  which  the  data  are  reported,  about  450  (^sec,  there  is  no 
influence  of  the  complicating  factors  due  to  unloading  and  the  prob- 
lem can  be  treated,  as  is  done  in  the  following  analysis,  as  a 
problem  of  loading  in  a semiinfinite  tube.  Material  properties 
like  shear  modulus  and  tangent  modulus  for  the  plastic  portion  of  the 
stress-strain  wave  were  taken  from  the  1967  report  of  the  authors. 

Value  of  material  density  for  test  specimens,  not  reported,  was  taken 
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from  tables  as  0.000871  lb  - sec  /in.  . 

2 . 2 Equations  of  Motion  and  Compatibility 

The  equation  of  motion  for  the  propagation  of  a torsional 
loading  wave  of  plastic  deformation  in  a thin-walled  tube  is  given 
as  follows  with  no  restriction  to  small  strains  implied. 


9t  Sv 

8x  = P "5t 


(2.1) 


The  compatibility  condition,  describing  the  continuity  of  tangential 
displacement  U(x,t)  of  any  material  point,  is  given  by 


_ 9v 
)t  ~ 


au.  . 


where  y (=  g^)  is  the  engineering  shear  strain. 


(2.2) 
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2. 3 Constitutive  Laws  Considered 

Three  constitutive  equations  will  be  examined  in  turn. 

(i>  I?  = I Hr  ♦ ltT  - f(v>]  • <2-3) 

This  equation,  incorporating  a linear  overstress,  was  used  by  several 
authors  as  noted  in  Chapter  1.  However,  Bianchi  (1963)  , employing  this 
same  equation  to  analyze  incremental  strain  data  in  longitudinal  impact, 
recommended  the  use  of  an  exponential  function  of  overstress  for  a 
better  fit.  This  is  tried  next. 

(ii)  It  = i if  + I • (2,4) 

(iii)  As  against  the  above  two  Malvern-type  semilinear  equations, 
we  chose,  finally,  a general  quasilinear  equation  of  the  type  used 
by  Cristescu  (1972). 

^=[l  + tp(T,Y)]|l+|[T-f<Y)]  (2.5, 

where  the  function  cp  has  to  obey  some  restrictions  mentioned  later. 


2 . 4 Interior  Differential  Equations 
and  Associated  Characteristics 

Equations  (2.1),  (2.2)  and  any  one  of  the  equations  (2.3), 

(2.4),  and  (2.5)  can  be  compactly  written  in  the  following  matrix 

form 


(2.6) 
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where  subscripts  t and  x denote  differentiations  with  respect  to  them 
and 

cp  = 0,  i|f  = — [T  - f (y)  ] for  equation  (2.3), 

cp  = 0,  t = | [exp  ( (t  - f(y)^4)  - for  equation  (2.4), 


cp  = cp(T,y),  + = - f (y) ] 


for  equation  (2.5) 


The  theory  for  systems  of  first-order  quasilinear  partial  differential 
equations,  as  represented  by  equation  (2.6),  is  given  in  Jeffrey  and 
Taniuti  (1964) . The  characteristic  wave  speeds  associated  with  such 


a system  are  the  roots  of  the  determinantal  equation 

pc  1 


Det 


c(9+i) 


= 0, 


(2.7) 


giving  c 


= (?+£)) 


with  the  characteristic  equation  given  by 


dx  = ± cdt  . (2.  8) 

(It  may  be  noted  that  for  equations  (2.3)  and  (2.4),  where  cp  = 0, 

c becomes  a constant,  equal  to  the  elastic  shear  wave  speed  «/G/p . ) 

T 

The  null  vector  i ( X , J?  ) for  equation  (2.6),  which  corresponds  to 

1 Z 

the  characteristic  wave  velocity  c,  is  given  by 


(1i  - V 


pc 


c(cp+i) 


= (0,  0) 


(2.9) 


Choosing  an  arbitrary  scale  factor  such  that  4 = 1 in  equation  (2.9), 


XT  = [-C(cp  + I),l]. 


we  get 


(2.10) 
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The  interior  differential  equation  corresponding  to  the  characteristic 
given  by  dx  = c dt  is  given  by 


which  is  rewritten  as 

dT  - p c dv  = - — dt , dx  = + c dt  . (2. 12) 

«P+±) 

For  the  negative  characteristic  given  by  dx  = -cdt,  we  replace 
(+c)  by  (-c)  in  equations (2. 12)  . 

2. 5 Boundary  and  Initial  Conditions 

Yew  and  Richardson  (1967)  report  incremental  shear  strain-time 
plots  for  three  strain  gauge  stations.  The  boundary  condition  for  the 
present  theoretical  solution  was  taken  to  be  the  incremental  strain-time 
data  for  the  first  strain  gauge  located  at  0.38  in.  from  the  impact  end. 

Initial  level  of  prestrain  developed  in  the  test  specimens  was 
not  mentioned  in  the  reports  (1967,1969).  Theoretical  solutions  com- 
puted with  two  different  levels  of  prestrain,  1 percent  and  7 percent, 
did  not,  however,  give  any  significant  difference  in  the  corresponding 
results.  Evidently  this  is  due  to  the  linearity  in  the  quasistatic 
stress-strain  relation  in  the  plastic  range  indicated  by  the  experiments 
(see,  for  instance,  the  1967  report)  and  used  in  the  theory.  The 
theoretical  solutions  presented  later  make  use  of  an  initial  strain  of 
1 percent  and  the  corresponding  stress  from  the  quasistatic  stress- 


strain  relation. 
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2 . 6 Linear  Overstress  - Numerical  Scheme 

Figure  2.1  a(i)  shows  the  solution  domain  in  the  x-t  plane. 

Mesh  point  types  involved  are  sketched  in  the  same  figures  with  labels 
(ii)  , (iii) , (iv) , and  (v) . The  unknowns  at  a typical  point  are 
stress  and  velocity  for  the  boundary  and  stress,  strain  and  velocity 
for  a general  interior  point.  These  are  determined  by  solving  the 
finite  difference  representations  of  equations  (2.12)  and  (2.3)  given 
as  follows  where  a and  b are  taken  for  the  equation  for  the  quasistatic 
stress-strain  relation,  = ay  + b. 

General  boundary  point  (Ref.  Fig.  2.1  a(iii)) 

Solution  is  known  at  points  B and  C.  Solution  at  point  T is 
obtained  from 

ST 'SC  + PC(VT  “V  = " AT  * K [0-5(ST  + sc>  "°.5a(GT+  Gc)  -b]  (2.13) 

St"SB_R(GT_GB)  = “ AT  * K [°-5(St+  Sb)  -°.5a(GT+  ty  -b]  (2.14) 

where  R stands  for  the  rigidity  modulus. 

Mesh  point  type  2.1  a(ii) , representing  the  first  point  to  be 
computed,  is  governed  by  equations  identical  to  those  above,  except 
that  AT  in  equation  (2.13)  is  replaced  by  0.5AT.  State  variables 
indicated  with  a subscript  C were  obtained  from  the  initial  condition 
for  this  case  of  no  step  function  boundary  loading. 

General  interior  point  (Ref.  Fig.  2.1  a(v)) 

* 

Solution  at  points  A,  B,  C are  known.  Solution  at  point  T is 


obtained  from 
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Figure  2.1  Schematics  of  solution  domains  and  mesh  point  types 
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i 

ST_SA"PC(VT  “V  = - AT  • K[0.5(ST+ SA)  -0.5a(GT+ GA)  -b]  (2.15) 

ST-Sc  + pc(VT-Vc)  = - AT  •K[0.5(ST+SC)  -0.5a(GT+Gc)  -b]  (2.16) 

St-Sb-R(Gt-Gb)  = - AT  • K[0.5(ST+  SB)  -0.5a(GT  + GB)  -b]  . (2.17) 

Mesh  point  behind  leading  wave  (Ref.  Fig.  2.1  a(iv)) 

The  difference  equations  are  the  same  as  equations  (2.15), 
(2.16),  and  (2.17)  except  that  AT  in  equation  (2.16)  is  replaced  by 
0.5AT.  G , S , and  V were  obtained  from  initial  conditions  in  this 

c c c 

case,  while  for  the  case  of  a leading  plastic  shock  wave,  evaluation 
of  an  integral  is  necessary  (Malvern,  1951). 

The  computations  were  done  with  a constant  time  step  of 
AT  = 1.42775  )i,sec.  This  time  step  was  chosen  so  that  the  third  strain 
gauge  fell  on  an  exact  mesh  point,  while  for  the  second  gauge  a Lagrange 
interpolation  was  necessary.  Initially,  boundary  data  were  read  at 
larger  intervals  and  a linear  interpolation  scheme  was  used  for  the 
intermediate  points.  The  data  so  generated  were  next  passed  through 
a standard  smoothing  routine.  Convergence  of  the  computations  was 
checked  by  seeing  that  halving  the  step  size  did  not  produce  any  observ- 
able change  in  the  solutions.  Results  are  discussed  in  Section  2.9. 

2.7  Exponential  Overstress  - 

An  Integro-Dif f erential  Approach 

Interior  differential  equations  (2.12)  for  the  exponential 

overstress  type  constitutive  law,  equation  (2.4),  can  be  explicitly 


written  as 
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dT  - pcdv  = - B dt  (2.18) 

where 

B = K |%xp  f |t  - f (y)j  /A ^ - lj  . (2- 19) 

The  finite  difference  form  of  this  and  equation  (2.4)  for  mesh  point 
types  given  in  Figure  2.1  a,  written  in  a predictor— corrector  format, 
results  in  three  nonlinear  algebraic  equations  in  the  case  of  an  inter- 
ior point.  Solutions  of  these  equations  computed  by  the  method  of 
Gauss-Seidel  iteration  failed  to  give  any  reasonable  match  with  the 
experimental  data.  The  theoretical  solutions  that  could  be  computed 
without  convergence  difficulties  always  showed  much  more  strain— rate 
effect  than  the  experimental  results.  Strain-rate  insensitivity  is 
meant  to  be  increased  by  increasing  K in  equation  (2. 19) , but  as 
Efron  (1964)  noted  earlier,  the  sufficiency  condition  for  convergence 
of  Gauss-Seidel  iterations,  places  a restriction  on  the  product  (K*AT) , 
such  that  the  higher  the  value  of  K,  the  smaller  is  the  step  size  AT 
that  permits  convergence. 

The  predictor-corrector  solution  by  Gauss-Seidel  iteration  was 
finally  abandoned  because,  even  for  the  smallest  step  size  that  could 
be  used  without  requiring  excessive  computer  time,  the  theoretical 
solution  still  showed  large  rate-dependence.  Apart  from  this  aspect, 
it  was  found  that  the  exponential  function  in  equation  (2.19)  was 
particularly  sensitive  to  a multiplicity  of  inflection  points  in  the 
boundary  data  that  remained  even  after  using  a smoothing  routine. 

It  was  thus  decided  to  use  a method  that  would  employ  integration 
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rather  than  differentiation.  This  would  take  care  of  any  disturbance 
at  the  boundary — integration  being  itself  a smoothing  operation — and 
as  later  turned  out  to  be  the  case,  would,  by  eliminating  one  of 
the  unknown  variables,  improve  the  picture  of  convergence  itself. 

This  is  the  integro-dif f erential  approach  discussed  below. 

Consider  the  following  ordinary  nonlinear  differential 
equation,  which  represents  the  constitutive  law  assumed  for  a mate- 
rial point. 

T - Gy  = - k j^exp  ((T-f(y))/A)  - lj  . (2.20) 

The  solution  to  the  above  equation  without  the  nonlinear  term  of  the 
right-hand  side  may  be  written  as 

T = G(y^y  ) + k(t-t  ) + p,  (2.21) 

where  p,  = const.  The  solution  to  the  complete  equation  (2.20)  may 
then  be  given  by  the  variation  of  parameter  technique  as  (such  an 
artifice  is  used,  for  example,  in  deriving  the  Krylov-Bogoliubov 
equations  of  nonlinear  oscillation  theory) 

T = G(y-yo)  + k(t-tQ)  + p,(t)  . (2.22) 

Substituting  equation  (2.22)  in  equation  (2.20),  we  get 

- p e“^/A  = k £exp  Q-  [G(y-YQ)  + k(t-tQ)  -f(y(t))])  J . (2.23) 

Recognizing  that  the  left-hand  side  in  equation  (2.23)  can  be 

written  as  A ~ (exp  (-u/A))  and  integrating  both  sides  between  t 
dt  o 

and  t,  we  get  [notice  that  p(&)  = T(0)  = T ], 

to  *0 


26 


M-(t) 


= - A ton  |_exp  (-  -2.)  + - / exp  ^ [G(y(s)  - Yq) 


+ k(s-tQ)  - f (y ( s ) ) ] J ds J . 


(2.24) 


Equation  (2.24)  substituted  in  equation  (2.22)  completes  the  solution. 
For  the  sake  of  brevity,  we  denote 


and 


1 r -] 

g(s)  = j jjKy(s)  - yQ)  + k(s-tQ)  - f (y(s))J 


i|f  = k(t-tQ)  - - A B/rX e To'/A  + ^ f eg^  ds^j 

*"0 


(2.25) 


Substituting  equation  (2.25)  in  equations  (2.24)  and  (2.22),  we  get 
an  explicit  representation  for  T in  terms  of  y, 


T — G(y  - y ) + T + llr  . 
‘ 1 o o 


(2.26) 


On  using  this  in  equation  (2.1),  the  equation  of  motion  is  now  written 
as  the  following  integro-diff erential  equation 


- G 


Bv  d[i|r] 

+ p Jt  = 


(2.27) 


where  i|f  is  given  in  terms  of  the  integral  in  equation  (2.25).  Equations 
(2.27)  and  (2.2)  now  constitute  a complete  set  of  functional  equations. 

In  these  equations,  if  we  treat  i|f  as  a dependent  variable  like  y and  v, 
we  observe  that  over  any  solution  domain  there  may  exist  lines  (character- 
istics) across  which  the  partial  derivatives  of  y and  v may  have  discon- 
tinuities, but  those  of  \Jr  (embodying  an  integral)  will  always  be  contin- 
uous. Thus,  while  looking  for  possible  discontinuity  lines  that  are 


\ 
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characteristics,  we  can  ignore  the  presence  of  the  term  involving  \|r  and 


7 

treat  equations  (2.2^)  and  (2.2)  simply  as  a system  of  first-order 


quasil inear  partial  differential  equations  as  before.  Thus  we  first 


write  them  in  matrix  form  as 


V J 


(2.28) 


The  characteristic  wave  speeds  are  obtained  from  (2.28), 


pc 


det 


= 0,  giving  c = 


G 


(2.29) 


1 c 


The  null  vector  (4  ,2  ) for  equation  (2.28)  which  corresponds  to  the 
characteristic  wave  velocity  c is  given  by 


(0,0). 


(2.30) 


Choosing  an  arbitrary  scale  factor  such  that  4 = 1,  we  get 


(2.31) 


The  interior  differential  equation  for  the  characteristic  given  by 
dx  = + c dt  is  then 
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which  simplifies  to 

dv  - c dy  = — dx  , dx  = + c dt  . (2. 33) 

T pc  ox 

For  dx  = - c dt  , (+c)  is  replaced  by  (-c)  in  equation  (2.33). 

Equation  (2.33)  is  integrated  next. 

2.7.1  Numerical  Scheme 

The  solution  domain  is  the  one  shown  in  Figure  2.1  b(i). 
Typical  mesh  points  are  shown  in  Figures  2.1  b(ii)  and  2.1  b(iii). 
Finite  difference  form  for  the  equation  (2.33)  is  written  as  follows: 


Boundary  point  (Ref.  Figure  2.1  b(ii)) 


V - V + c (G-G) 

= ±( 

(-AX) 

(2.34) 

T C T C 

where 

p C ' 

^ — ax  / 

i]r„  = k(t  -t  ) - T 
M Mo  o 

- A 2m 

1 exp  ^ 

(— T / A)  J + -T-  f exp  (g(s)) 
o / A 

ds 

o 


(2.35) 

Assuming  that  the  state  variables  V,  G,  and  (for  the  sake  of  explic- 
itness) \|f  are  known  at  C,  V can  be  determined  readily  from  above. 

is  also  directly  evaluated  from  equation  (2.26).  Next  we  take 
up  a general  interior  point  (Figure  2.1  b(iii))  where  the  solution 


is  known  at  A,  B,  C. 


29 


General  interior  point  (Ref.  Figure  2.1  b(iii)) 


V,  - V - c(G  -G,) 
T A T A 


(2.36) 


Vm  - V„  + c ( G -G_)  = 
T C T C 


_1_ 

pc 


L^J  (-4X> 


(2.37) 


Adding  the  two,  we  get 

vT  = °-5rvA  + vc  + c(Vga)  + (^c_V/pc]  • (2,38) 

Subtracting  equation  (2.36)  from  equation  (2.37)  we  have 

Gt  = 0.5[Gc  + Ga+  (Vc-Va)/c  + (fc+tA-21rM)/pc2]  . (2.39) 

In  equation  (2.39),  \jr  involves  (see  equation  (2.35))  the  unknown 

M 

strain  at  the  point  M for  which  we  iterate  as  follows. 


= k(t  + 0.  5 At  -t  ) - T - A 2n  \ exp  f(-T  /A)') 

M B oo  L\o/ 

tfi 

+ ^ | J exp  (g ( s ) ) ds  + 0 . 5 A t exp  (g(tfi))j-  J • (2.40) 
*o 

Substituting  equation  (2.40)  in  equation  (2.39),  we  get  a first 
approximation  for  G . Denoting  this  by  and  noting  that  the 

first-order  strains  at  M and  T may  be  related  to  the  strain  at  B as 


Y 


(1) 

M 


= 0.5 


Yn 


(1) 


+ y 


B 


we  have 
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k(tM'to)  " To  “ A^[exp  ((_To/A))  + | { / exp  (g(s))  ds 

o 

+ 0. 25A t (^exp  (g(tB>)  + exp  j [G(y^-  yQ) 

♦ *<vv-fO))]-  <2-41) 

Equation  (2.41)  is  now  substituted  in  equation  (2.39)  to  get  G . 

The  process  is  repeated  till  the  difference  between  two  successive 
iterations  on  G^  is  within  the  tolerance  limit.  Knowing  G^,  we  eval- 
uate i|r  from  equation  (2.25)  and  then  equation  (2.26)  gives  directly. 

The  superior  stability  property  of  this  scheme  of  computation 
over  that  of  the  Gauss-Seidel  method  can  be  readily  appreciated  from 

9 

the  fact  that  while  the  values  of  k = 10  and  A = 250  in  equation  (2.4) 
could  be  used  in  the  computations  by  the  integro-dif f erential  approach 
for  a step  size  of  AT  = 1.1422  jisec,  the  Gauss-Seidel  technique  could 

g 

admit  parameter  values  only  up  to  k = 3 x 10  and  A = 1000  for  a step 
size  of  AT  =’0.5711  (j£ec.  The  only  problem  with  the  integro-dif f eren- 
tial scheme  is  that  the  integrand  in  equations  (2.35)  and  (2.41)  soon 
exceeds  the  largest  number  that  the  machine  can  handle.  This  diffi- 
culty was  overcome  by  developing  a recursive  procedure  for  evaluating 
the  integrals  in  multiple  time  spans,  as  indicated  in  Appendix  A where 
a scheme  is  also  given  of  the  actual  computer  program  used  for  the 
integro-diff erential  approach.  The  machine  time  required  for  this 

problem  was  2.7  minutes  in  an  IBM  360/65  model  for  computations  up  to 

% 


370  psec. 
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2.8  Quasilinear  Model 

We  now  consider  the  behavior  of  the  theoretical  solution  for 
the  quasilinear  constitutive  model, 


2.8.1  Remarks  Concerning  the  Choice 
of  the  Function  c p 

The  function  cp  in  equation  (2.5)  is  associated  with  the 
concept  of  instantaneous  plasticity,  an  assumption  that  has  no  basis 
in  metal  physics.  Thus  the  function  cp  was  introduced  in  the  liter- 
ature [see  Cristescu  (1967)]  on  purely  phenomenological  grounds,  in 
order  to  give  a better  description  of  the  dynamic-plastic  response  of 
a material.  As  such,  the  procedure  for  choosing  the  form  of  the  func- 
tion cp  in  the  above  equation  is  rather  empirical.  We  will  follow  a 
scheme  that  has  been  outlined  by  Cristescu  (1972).  In  this  method, 
one  chooses  a simple  function  as  a first  approximation  to  the 
"instantaneous  response  curve"  (Cristescu,  1967)  where  the  slope  is 
taken  to  be  l/(cp+l/G),  and  then  by  adjusting  on  the  basic  form  for  cp 
a good  representation  of  the  rising  part  of  the  strain-time  curve  up 
to  an  inflection  point  can  be  obtained. 

The  basic  form  chosen  for  the  instantaneous  response  curve  in 
this  case  was 


(2.5) 


with 


Cp  = 0 if  (i)  T < f (y) 


or  (ii)  ^ < 0 . 


T = A J Y + Y1 


(a) 
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where  A and  are  the  parameters  to  be  determined.  The  desired 
instantaneous  response  curve  consistent  with  equation  (2.5)  is  given 
by  (obtained  by  setting  dt  -*  0 in  equation  (2.5)) 


dy  = (9  + p) 


dT 


(b) 


Calculating 
we  have 


from  both  equations  (a)  and  (b)  and  equating  them, 
dy 


(c) 


Now  setting  the  conditions  that  the  instantaneous  curve  passes 
through  the  initial  stress-strain  point  and  that  it  has  an  initial 
slope  equal  to  the  elastic  slope,  we  have 


T = A Jy  + y (d) 

o 'o  T1 

A , . 

——————  = G . ( e) 

2 ^y0+  y1 

Solving  equations  (d)  and  (e)  for  A and  y^  and  substituting  in 
equation  (c) , we  get 


2 J y -y  + T /2G 
' 'o  o 

o 


(f) 


This  is  the  first  approximation  form  to  the  instantaneous  response  mod- 
ulus. Several  trials  were  needed  to  arrive  at  a suitable  modification 
of  the  above  form  before  a reasonably  good  fit  was  obtained  with  the 
experimental  strain-time  data.  In  concurrence  with  Cristescu  (1972), 
it  was  observed  that  a higher  value  for  the  function  cp  had  the  effect 
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of  decreasing  the  length  of  the  (Ay  - 1)  curve  in  the  theoretical  solu- 
tion beyond  its  inflection  point.  The  form  of  the  coefficient  func- 
tion for  which  results  are  reported  in  this  thesis  was  chosen  with 
the  following  modification  on  equation  (f)  above. 


cp  = 


2 */u( y - yQ)  + Tq/2G 
*/2GT 


1 

G 


(2.42) 


with  a = 1500. 

• A final  restriction,  proposed  by  Cristescu  (1972) , on  the  form 
of  cp  in  equation  (2.5)  is  that  the  slope  of  the  dynamic  curve  must  lie 
between  those  of  the  instantaneous  curve  and  the  relaxation  boundary, 
such  that  unloading  does  not  occur  with  increasing  stress.  For  the 
present  case  then,  with  the  relaxation  boundary  given  by  T = f(y), 


f/(Y>  <5pTT<G- 

A predictor-corrector  solution  of  equation  (2.5)  with  cp  given  by 
equation  (2.42)  showed  no  crossing  of  the  relaxation  boundary  for  the 
increasing  strain-time  data  at  the  boundary,  and  thus  the  above 
inequality  condition  was  not  violated. 

It  may  be  noted  in  passing  that  Myers  and  Eisenberg  (1973)  have 
recently  developed  a theoretical  model  for  constructing  the  function  cp 
under  combined  stress  conditions.  Appropriate  reductions  from  these 
general  quasilinear  equations  yield  all  the  special  forms  for  cp  used 


by  Cristescu  (1972). 


* 
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2.8.2  Numerical  Solution  by  Second- 

Order  C.I.R.  Scheme 

For  the  quasilinear  model  we  now  consider  integrating 
equation  (2.5) — with  equation  (2.42)  incorporated  in  it — together 
with  the  interior  differential  equation  (2.12).  However,  in  this  case 
the  associated  characteristics,  though  still  given  by  equation  (2.8) 
in  form,  are  in  effect  functions  of  the  unknown  strain  at  any  point, 
and  thus  cannot  be  known  in  advance  as  before.  A second-order  accur- 
ate difference  scheme  for  the  numerical  solution  of  such  systems  of 
quasilinear  first-order  hyperbolic  partial  differential  equations  has 
recently  been  developed  by  Ranganath  and  Clifton  (1971),  and  called 
(by  them)  a second-order  Courant,  Isaacson,  Rees  method.  In  essence, 
the  method  consists  of  a two-step  integrating  routine.  First,  a pro- 
visional solution  is  obtained  by  the  standard  (first-order)  Courant, 
Isaacson,  Rees  procedure  (see  Jeffrey  and  Taniuti,  1964).  Then  this 
solution  is  used  in  evaluating  the  average  value  of  the  coefficients 
in  the  interior  differential  equation,  while  the  interpolation  for  the 
mesh  variables  at  the  intersection  of  the  characteristic  with  the  base 
line  of  constant  time  is  done  by  Taylor  expansions,  keeping  up  to  the 
second  derivative.  We  will  use  this  method  in  the  present  case. 

The  solution  domain  is  schematically  represented  in  Figure 
2.1  c(i).  The  computations  were  done  for  a constant  time  for  any  row 
of  points.  The  basic  rectangular  grid  was  constructed  for  any  time 
step  At  such  that  AX  = CAt  (C  = elastic  speed),  and  this  guaranteed 
stability  of  the  computations.  There  were  five  types  of  mesh  config- 
urations involved,  as  shown  in  Figures  2.1  c(ii)  - (vi) . Of  these, 


35 


only  a general  interior  point  computation  employs  the  complete  algo- 
rithm as  proposed  by  Ranganath  and  Clifton  (1971) . Points  at  the 
boundary  or  those  lying  just  behind  the  leading  wave  front  needed  a 
slight  modification  because  of  the  different  mesh  geometries  involved. 
The  first  three  computed  points,  as  indicated  in  Figures  2.1  c(ii)  - 
(iv) , are  treated  by  conventional  iteration  schemes.  We  will  outline 
here  only  the  computational  schemes  for  two  basic  mesh  types  in  any 
general  row,  while  the  flow  chart  for  the  computer  program  used  is 
given  in  the  Appendix. 


General  boundary  point  (Ref.  Figure  2.1  c(v)) 

In  this  and  the  subsequent  mesh  types  discussed,  we  rewrite 
equation  (2.12)  as  follows, 


obtaining 


dr  + p c dv  = - p c i|r  dt 


dx  = ± c dt 


(2.43) 


c = ——————  . (2.44) 

Vp  (9+1/ G) 

In  the  computations  mentioned  below,  cp  is  set  equal  to  zero 
(i)  whenever  T at  an  unknown  point  is  computed  in  the  course  of  iter- 
ations as  less  than  the  quasistatic  stress  for  the  corresponding  strain, 
or  (ii)  if  the  computed  stress  is  less  than  the  stress  in  the  previous 
time  step  for  any  x. 

Superscripts  (1)  and  (2)  refer  to  first-  and  second-order 
solution,  respectively,  in  the  following.  Subscripts  refer  to  Fig- 
ure 2.1  c(v).  G is  known  from  the  boundary  condition. 
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- GD  = 0.5^  + CpD)  (S^2)-  SJ  +0.5AT[|(ST-aGT-b)  + BD  | (2.45) 


T B 


T TB  T B 


J 


where  9 = (cp  + 1/G)  here  and  elsewhere,  R = rigidity  modulus, 
BB  = |(SB"aGB'b)- 


c(G  ) 

6 = G + (G  - G ) , 

R B C C B ’ 


C - elastic  shear  wave  speed 


X = 1 - 


c(Gt)  + c(fiR) 
2~C 


gr  = Gc  - °-5Ugd-gb>  + °-5X  <gd-2gc  + V 


Sr  = Sc  - 0.5USd-Sb)  + 0.5  X (Sd-2Sc+Sb) 


Vr=Vc-°.5X(Vd-Vb)  + 0.5  X (Vd-2Vc  + Vb) 


C = 0. 5 C(G„>  + C(G„) 
m V T R 


Si2,-SR  + P<v<2)-VR,. 


- 0.  5C2  AT  (b(2)  + B 
m \ T R 


(2.46) 


B^2^  =^(s^,2)-aG  - b');  B = -(s  - aG  - b 
T R\  T T /’  R R\R  R 


(2)  (2) 

Equations  (2.45)  and  (2.46)  give  S^,  , directly. 


General  interior  point  (Ref.  Figure  2.1  c(vi)) 

The  first-order  solution  is  constructed  as  follows: 


c(gb} 

SL=SB-—  (SB-V 


C(V 

SR  = SB  + C (SC  _ V 


C(G) 
v B 

V = V — (V  -V  ) 

L B C B A 
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V_  = 


C(V 

v„  + — - — (v  - V ) 


R B C C B 


C(GB) 

GL=  Gb'  — ^B-V 


C(G  ) 

v B 

G = G + (G  - G ) 

R B C V C B 


cii  v r c i v ~ I 2 

ST  " SL  -PC(V  LVT  -VJ=-PC  (V  BBAT 


(2.47) 


S^1)  - Sc  + pc(GB)  [vrJ1)  -vj  = - pc2(Gn)  B0  AT 


R 


B B 


(2.48) 


G(1)  - G = cp  (S(1)  -S  ) + B At. 
T B yB  T B B 


(2.49) 


Solution  of  equations  (2.47),  (2.48)  and  (2.49)  gives  the 
values  G^,^ . Second-order  solution  is  computed  as  follows: 


\ 2c{_ 


c(G^,1))  + c(Gl) 


] 


2C 


(1)  ^ 
c(G^±;)  + c(Gr) 


J 


G,  = 


- gb-°-5Vgc-ga)  + °-!W20b+0a) 


VL  * VB  -°- 51<VC-VA>  + °-5V(VC-2VB+V 


Bl=  Sb-0.5X(Sc-Sa)  ♦ 0.5Xl(Sc-2Sb+Sa) 


G„  = 


Gb+  °.5XWc-Ga)  + 0.5Xr(Gc-2Gb+Ga) 


\ = VB  + o^To-V  + °-5Xr<vc-2vb  + va) 
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SE  = SB  + °-5HSC-SA>  + °-5XR<Sc-2SB+SA) 


mL 


= 0.5  [c(G«)  ♦ coy}  = 0.  5 [c(G^1)  ♦ c(Gr)] 


, (2) 


- S_  - pc  t(v!2)  -V)  = -pc2  (0.5  AT)  jB  +B  1 (2.50) 

L mL  T L mL  L T(l)  LJ 

BL  = ! [\  - aGL  - b]  ; ' B„<1)  = I [4"  - *4”  - bj 


s-2>  -s. + p<=~»  W2>  -O  = -pcl(o-5iT>  [BT(D  +BJ 


R mR  \ T R 


(2.51) 


B„  = 


S_  - aG_ 


R R L R R 


-] 


(2)  — - (2) 

GT  " GB  = °-5%  + CP  (1))(St  "V  + °-5AT  (BB  + B (1)}  (2,52) 

where  cp  = cp(G^^).  Solution  to  the  last  three  equations  gives 

s<2),  g<2),  v<2\ 


Mesh  Point  Lying  Just  Behind 
the  Leading  Wave 

In  the  problem,  all  points  on  the  leading  wave  are  in  the 

initial  rest  state  of  non-zero  plastic  shear  stress  T . The  computa- 

o 

tion  of  mesh  points  of  the  type  shown  in  Figure  2.1  c(iv)  is  done  on 
this  basis.  The  iterations  cannot  follow  the  interpolation  pattern 
outlined  above  because  of  the  peculiar  geometry,  and  instead  succes- 
sive substitutions  are  done,  with  linear  interpolations,  to  obtain  the 
solutions  as  in  the  case  of  the  mesh  point  shown  in  Figure  2.1  c(iii). 
The  first  iteration  in  such  cases  is  always  done  by  the  standard  first- 
order  C.I.R.  method. 
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All  the  above  three  mesh  types  involve  no  iteration  beyond  two, 
as  per  the  scheme  proposed  by  Ranganath  and  Clifton  (1971);  however,  in 
the  actual  program,  use  was  made  of  greater  numbers  of  iterations  on 
the  same  lines,  in  case  the  desired  accuracy  (between  successive  iter- 
ations for  strain)  in  the  solutions  was  not  obtained.  For  mesh  points 
with  stress  states  very  close  to  the  relaxation  boundary,  the  stress 
computed  by  the  second-order  equations  was  often  found  to  be  less  than 
the  quasistatic  value.  In  such  cases  the  next  iteration  was  computed 
on  the  basis  of  elastic  equations.  If  the  stress  so  obtained  was  still 
less  than  the  quasistatic  stress , then  this  result  was  accepted;  other- 
wise, a "crossing  routine"  as  described  by  Cristescu  (1970)  was  used 
to  obtain  the  stress  for  the  point.  The  mesh  size  chosen  was  0.5711 
|isec  in  time  step.  Halving  of  this  step  size  did  not  produce  any  sig- 
nificant change  in  results.  The  total  computer  time  on  the  IBM  360/65 
machine  was  7.6  minutes  for  a loading  duration  of  370  psec. 

2. 9 Results  - Comparison  with  Experimental  Data 
and  Rate-Independent  Solutions 

Figures  2.2,  2.3,  and  2.4  present  the  incremental  strain-time 
plots  obtained  on  the  basis  of  the  three  constitutive  models  examined. 
Also  shown  in  these  figures  are  the  corresponding  experimental  results 
(Yew  and  Richardson,  1969)  and  solutions  based  on  a rate-independent 
theory  using  the  quasistatic  stress-strain  relation  as  the  dynamic 
curve. 

% 

Figure  2.2  corresponds  to  the  linear  overstress  type  semilinear 

0 

model,  equation  (2.3),  with  a value  of  k = 8.5  x 10  . It  is  seen  that 
though  the  third-gauge  response  is  well  described  from  about  0.20  per 


Rate-independent  solution 
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Figure  2.2-  Incremental  strain-time  plots  for  the  linear  overstress  model 
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cent  incremental  strain,  the  solution  at  the  second  gauge  consistently 
lags  the  experimental  data.  The  extremely  small  incremental  strains 
as  well  as  those  above  1.0  per  cent  are  not  well  described,  but  while 
very  low  strain  measurements  are  always  suspect,  the  large  strain 
measurements,  particularly  at  the  second  gauge,  seem  to  indicate  a 
delayed  response,  as  is  known  to  occur  with  strain  gauges  in 
detecting  large  plastic  strain.  In  any  case,  it  is  seen  that  the  rate- 
dependent  solution  gives  a better  overall  description  than  a rate- 
independent  solution  based  on  the  quasistatic  stress-strain  law. 

Figure  2.3  gives  the  incremental  strain  plots  for  the  expo- 
nential overstress  type  semilinear  model.  With  smaller  level  strains 
propagating  faster  and  the  higher  level  strains  slower  than  in  the 
experiments,  the  solution  presented  is  indicative  of  a stronger  rate- 
effect  than  that  represented  by  the  measurements.  The  situation  could 
be  improved  with  a higher  value  of  k in  equation  (2.6),  but  even  with 

9 

the  present  one,  k = 10  with  A = 250,  the  integral  in  equation  (2.25) 
had  to  be  evaluated  in  13  parts  to  prevent  overflow,  and  further 
increase  in  k would  have  made  the  program  too  unwieldy.  Even  so,  the 
general  agreement  between  the  present  rate-dependent  solution  and  the 
experimental  measurement  is  greater  than  in  the  case  of  the  rate- 
independent  solution. 

Figure  2.4  shows  the  results  for  the  quasilinear  model.  Data 
at  the  second  gauge  are  matched  fairly  well,  with  incremental  strains 
above  .012  lying  between  the  experimental  values  and  the  rate- 
independent  solution.  Data  at  the  third  gauge  perhaps  agree  as 
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Figure  2.3  Incremental  strain-time  plots  for  the  exponential  overstress  model 
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Figure  2.4  Incremental  strain-time  plots  for  the  quasilinear  model 
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as  well  as  that  with  the  exponential  overstress  model.  Again  the  gen- 
eral agreement  between  the  theory  and  the  experiment  is  better  than 
with  the  rate-independent  solution.  Also,  among  the  three  theoretical 
solutions,  this  perhaps  gives  the  best  overall  fit  to  the  data. 

Figure  2.5  shows  a stress-time  plot  at  the  boundary  for  the 
three  constitutive  models.  It  may  be  noted  that  the  boundary  stress 
could  not  well  be  obtained  from  these  experiments,  at  least  not  with 
our  choice  of  boundary  at  the  first  gauge.  It  is  seen  that  with  respect 
to  the  rate-independent  solution,  that  is,  the  quasistatic  stress  val- 
ues, all  the  three  solutions  show  a slightly  higher  stress  initially 
due  to  the  so-called  overstress  in  the  rate-dependent  theory.  When 
the  strain  rate  becomes  very  small  at  the  boundary,  all  the  rate- 
dependent  solutions  give  the  same  result  as  the  rate-independent 
theory  based  on  the  quasistatic  curve. 

2.10  Comparative  Study  of  the  Three  Models 
and  Approaches 

The  rate-independent  solutions  plotted  in  Figures  2.2,  2.3,  and 
2.4  for  the  second  and  third  strain  gauge  stations  show  that  the  smaller 
level  strain  propagated  faster  in  the  experiment  than  such  a theory 
predicted,  while  for  the  higher  level  strains  the  situation  was  reversed. 
This  was  established  by  Malvern  (1951)  to  be  indicative  of  the  strain- 
rate  effect.  Thus  it  is  no  surprise  that  all  of  the  three  models  con- 
sidered showed  better  agreement  with  the  experimental  results  than 
a rate-independent  solution.  The  difference  between  the  three  solutions 
does  not  seem  to  be  large,  though  the  quasilinear  model  seems  to  give 
a slightly  better  overall  fit  to  the  data  available. 


11441 
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Figure  2. 5 Stress-time  plot  at  the  boundary 
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Regarding  the  methods  applied  for  solution,  the  linear  over- 
stress model  gives  the  least  trouble  in  computation.  Because  the 
quasistatic  stress-strain  law  used  was  linear,  the  predictor-corrector 
equations  could  be  directly  solved,  and  unless  a very  high  value  of  k 
was  used  in  equation  (2.3),  no  numerical  oscillation  occurred  in  the 
computations  if  the  step  size  was  not  too  large. 

The  integro-dif f erential  approach  definitely  proved  itself 
to  be  a powerful  technique.  The  method  can  handle  even  rough  boundary- 
data  and  works  successfully  over  parameter  ranges  for  which  iteration 
schemes  like  the  Gauss-Seidel  method  fail  due  to  convergence  restrictions. 

The  second-order  Courant,  Isaacson,  Rees  (C.  I.R.  ) method,  as 
proposed  by  Ranganath  and  Clifton  (1971)  requires  only  two-step  inte- 
gration with  no  iterations.  For  the  step  size  used  in  the  present 
problem,  0.5711  (jsec,  which  is  not  large  considering  the  duration  of 
impact  for  which  the  solution  was  sought,  such  a procedure  led  to  numer- 
ical instability.  The  situation  was  corrected  by  making  a provision 
for  iterations  up  to  five  in  number,  starting  with  the  second-order 
solution.  The  method  in  general  required  the  maximum  computer  time 
among  the  three  solutions. 

2.11  Conclusions 

1.  For  the  material  used  in  the  wave  propagation  experiment 
under  consideration,  a rate-dependent  theory,  with  appropriate  values 
for  the  associated  constants,  definitely  ^ives  a better  agreement  with 
the  experiment  than  a rate-independent  theory  based  on  the  quasistatic 


stress-strain  relation. 
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2.  Different  rate-type  laws,  or  the  same  law  with  different 
values  for  its  constants,  can  apparently  give  the  same  incremental 
strain-time  solution.  This  means  that  such  measurements  alone  are 
no  definitive  indicators  of  the  constitutive  behavior.  Additional 
kinds  of  data  are  needed  to  determine  which  constitutive  formulation 
is  most  applicable  in  a variety  of  dynamic  loading  situations. 

3.  For  the  numerical  solution  of  nonlinear  stress  wave  prob- 
lems, both  the  integro-diff erential  approach  and  the  second-order  C.  I.R. 
method  represent  powerful  techniques  that  are  superior  to  the  conven- 
tional iteration  schemes  like  the  Gauss-Seidel  method,  in  that  they 
work  successfully  over  a parameter  range  for  which  conventional  iter- 
ations would  not  converge.  This  speaks  of  the  superior  stability 
property  of  these  methods. 


CHAPTER  3 


INCREMENTAL  LOADING  WAVES  OF  COMBINED  STRESS 

In  this  chapter  we  consider  a problem  on  combined  longitudinal 
and  torsional  plastic  waves  in  a prestressed  tube.  A rate  dependent 
theory  is  formulated  based  on  the  assumption  of  a viscoplastic  poten- 
tial function  and  specific  hardening  mechanisms.  Numerical  solutions 
are  compared  with  the  experimental  data  and  rate-independent  solution 
of  Lipkin  and  Clifton  (1970). 

3. 1 The  Experiments  of  Lipkin  and  Clifton 

Lipkin  and  Clifton  (1969,1970)  made  an  experimental  study  of 
plastic  waves  of  combined  stresses  caused  by  a longitudinal  impact  of 
a pretorqued  tube.  A relatively  rate-insensitive  material,  annealed 
3003H14  aluminum,  was  chosen  for  the  specimen  which  was  a tube  with 
0.5  inch  (nominal)  outside  diameter,  0.049  inch  wall  thickness  and 
41  inch  length.  This  was  glued  to  a solid  elastic  transmitter  bar  of 
titanium  and,  just  before  impact,  the  transmitter-specimen  assembly  was 
subjected  to  a static  torque  sufficient  to  cause  plastic  deformation 
in  the  specimen  tube.  The  far  end  of  the  specimen  is  restrained  from 
rotation  by  a rear  mounting  assembly,  but  is  free  to  move  axially. 

A short  titanium  hitter  bar,  released  by  a gas  gun,  causes  the  impact 
on  the  free  end  of  the  transmitter  bar.  Strain  gauges  were  used  to 
measure  strain  at  three  locations  on  the  specimen  and  two  on  the 
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transmitter  bar.  Direct  measurement  of  the  velocity  of  the  hitter  bar 
was  done  by  means  of  an  electronic  counter,  started  and  stopped  as  the 
hitter  passed  two  accurately  spaced  photodiodes.  Particle  velocities 
at  the  transmitter-specimen  interface,  obtained  from  this  measurement 
as  well  as  the  elastic  computations  in  the  transmitter  bar,  showed  that 
excepting  for  a small  rise  time  these  were  essentially  constant  up  to 
a certain  time,  following  which  unloading  occurred.  Experimental 
results  reported  by  these  authors  (1969,1970)  were  up  to  a time  before 
which  no  unloading  wave  could  be  reflected  from  the  free  end,  and  so 
in  the  following  the  problem  is  analyzed  as  purely  of  loading. 

3. 2 Equations  of  Motion  and  Compatibility 
under  Combined  Stresses 

We  consider  the  following  material  coordinate  system  to  describe 
the  motion  due  to  impact  for  the  thin-walled  tube:  axis  X is  longi- 
tudinal, is  tangential,  and  X^  is  radial.  The  equation  of  motion, 

then,  under  the  assumption  that  stresses  are  independent  of  X and  X„ , 
can  be  written  in  the  reference  state  of  undeformed  configuration  as 

dT^  5v^ 

“5x^  = po  ~W  (a) 

O 

where  T is  a component  of  the  first  Piola-Kirchoff  stress  tensor, 

p is  the  density  in  the  reference  state,  and  v are  the  particle 
o i 

velocity  components  at  time  t.  In  the  present  study  we  neglect  T° 

13 

and  the  radial  inertia  p (dv„/3t) . The  relation  between  the  Piola- 

o 3 

Kirchoff  stress  tensor  and  the  more  common  Cauchy  stress  tensor  (see, 
for  example,  Malvern  (1969),  p.  223)  in  this  case  reduces  to 
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aH  " P_ 
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"dx 

1 m° 

3x^  T11 


dx 

1 T° 

ax.  12 


(b) 


J 


where  x^  = , x^ , U defining  the  current  spatial  coordinate 

(with  respect  to  the  reference  coordinate  system  at  undeformed  state) 
and  the  particle  displacement,  respectively,  and  p is  the  current 
density.  Equations  (b)  are  rewritten  as 


'N 


11 


CT.  „ = 


12 


f (i  + eii>  Tii 

o 


t (1  + 6ll)  T°2j 


(C) 


where  = Bu^/dX  . 

The  density  ratio  in  equation  (c)  can  be  written  from  the 
material  form  of  the  continuity  condition  (see,  for  example,  Malvern 
(1969),  p.  208)  which  reduces  in  this  case  to  (see  Goel  (1969)) 


— = U + e11)(l+e22)(l  + e33) 


(d) 


where  and  e33  are  the  non-zero  normal  strains  in  the  tangential 
and  radial  directions,  developed  due  to  the  Poisson  effect.  Substi- 
tuting equation  (d)  in  equations  (c) , we  have 


T11  " (1 + e22) (1 + e33)  Q11 


T12  (1  + e22) (1  + e33)  °12  • 


(e) 
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Replacing  for  T°  , T°2  from  equation  (e)  in  equation  (a), 


we  get 


^ [(1+  S22)(1+ e33}  CTlJ=  Po  “SF 

Cl2]=  Po 


9 r 


(1+e22)(1+e33} 


^2 
~5T  • 


(f) 


(g) 


Equations  (f)  and  (g)  represent  the  equations  of  motion  for  the  present 
problem  for  finite  strain.  If  we  restrict  ourselves  to  small  displace- 


ments and  strains,  then  e 


22 


°’  e33  ~ °'  X1  ~ xi>  po  ^P 


Writing 


CT11  aS  °12  as  T » as  u and  v2  as  v>  and  X1  as  x»  we  then  obtain 
the  equations  of  motion  for  infinitesimal  strains, 


So 

Su 

P St 

(3.1) 

St 

Sv 

^ = 

p St  • 

(3.2) 

This  is  the  equation  of  motion  that  we  will  use  in  the  subsequent 
analysis. 

The  compatibility  conditions,  arising  out  of  the  continuity  of 
the  functions  describing  the  displacements  (in  the  radial  and  tangen- 
tial directions)  with  respect  to  space  and  time,  are  given  by 


Su 

Se 

3S  = 

St 

(3.3) 

Sv  Sv 

"St 

(3.4) 

Su 


where  y 
in  the  X, 


(=  sr) 1 


is  the  engineering  shear  strain  for  displacement  Ur 


direction,  and 


e = 


dvi 
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3.3  Constitutive  Laws  Considered 

In  this  section  we  construct  constitutive  models  to  describe 
rate-dependent  plastic  behavior  under  combined  stress.  The  material 
is  being  assumed  to  be  elastic-viscoplastic  with  possible  work- 
hardening. 


3.3.1  Viscoplastic  Potential  Function 
and  Associated  Flow  Rule 

Rice  (1970)  introduced  the  idea  of  a viscoplastic  potential 
function  D under  the  assumption  that, for  any  slipped  state  of  disloca- 
tions, the  rate  of  slipping  is  governed  by  the  resolved  shear  stress  or 
by  the  local  forces  on  dislocation  lines.  A specific  form  satisfying 
path  dependence  and  implying  isotropic  and  kinematic  hardening  was 
given  as 

fl  = f [J*(£  - J3)  , k]  (a) 


where  ct  is  the  stress  tensor,  3 is  the  matrix  representing  the  coor- 
dinates of  the  center  of  the  yield  hypersurface  in  nine-dimensional 
* 

stress  space,  J is  the  second  invariant  of  the  tensor  represented  by 
(ct  - 3) , and  k is  a work  hardening  parameter.  With  the  assumption  of 
a viscoplastic  potential  function,  the  structure  of  the  associated  flow 
rule  for  the  case  of  only  two  non-zero  stress  components  (<r,T)  is 
[see  Rice  (1970)]  given  as 


"3~t  = ^ 

Ot  OT 


» 


(b) 


(c) 
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With  the  object  of  particularizing  these  equations  to  our  problem,  we 
assume  a viscoplastic  potential  function  of  the  following  form. 


where  in  this  case  of  only  two  non-zero  stress  components  (c,T) , 

= a*)2  + (r  - T*)2  (e) 

* * 

with  a , T indicating  the  coordinates  of  the  center  of  the  yield 
surface.  The  motivation  behind  assuming  this  particular  form  of  the 
viscoplastic  potential  function,  equation  (d) , is  to  match  the  asso- 
ciated flow  rules  with  those  derived  by  Perzyna  (1962)  from  a general- 
ization of  the  Hohenemser-Prager  equations.  Thus,  substituting 
equations  (d)  and  (e)  in  equations  (b)  and  (c) , we  have 


~5F  = 


2(ct-o*) 


exp 


(3.5) 


These  are  precisely  the  equations  that  are  obtained  from  Perzyna’ s 
constitutive  equations  for  rate-sensitive  plastic  materials  if, 

(1)  one  chooses  in  Perzyna' s equations  (see  Perzyna,  April  1962)  the 
functional  form  given  as  follows, 
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t-rn 


exp  a ^ — — 


V 


and  (2)  combined  translation  and  expansion  of  the  yield  surface  is 
understood.  To  get  the  Malvern-type  linear  overstress  equations 
corresponding  to  combined  stress,  one  would  consider  the  form 


Equations  (3.5)  and  (3.6)  will  be  used  in  the  following  as  the  govern- 
ing plastic  strain-rate  equations. 

3.3.2  Generalized  Combined  Hardening 
Formulation 

Rice's  theory  of  a viscoplastic  potential  envisages  combined 
isotropic  and  kinematic  hardening.  A method  of  apportioning  the  total 
hardening  into  the  isotropic  and  kinematic  modes  was  formulated  by  Goel 
and  Malvern  (1970)  for  their  studies  on  combined  stress  wave  propagation 
in  materials  exhibiting  nonlinear  stress-strain  law.  The  following 
development,  based  on  the  latter  work,  th&s  fits  in  the  scheme  of 


corresponding  to  which  equations  (b)  and  (c)  assume  the  form 


2 


Rice's  proposition. 
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Figure  3.1  shows  schematically  in  two  dimensions  the  processes 
that  are  implied  by  combined  hardening.  Ellipse  A denotes  the  orig- 
inal yield  surface  in  the  o—T  plane,  and  ellipse  B is  the  translated 
and  expanded  yield  surface,  a is  the  vector  in  the  ct-t  plane  repre- 
senting the  tensor  components  in  nine-dimensional  stress  space,  3 is 
similarly  the  vector  describing  the  position  of  the  translated  yield 
surface,  and  is  the  stress  tensor  referred  to  the  center  of  this 
translated  surface. 

In  analogy  with  isotropic  hardening,  the  equivalent  "effective 
stress”  to  cause  yield  in  uniaxial  tests  is  given  by  the  scalar  R (see 
Figure  3.1)  defined  as 


where  §*=§..-i§  6 

ij  ij  3 ^kk  ij 


and  %..  = a - B 

ij  iJ  ij 

When  the  hardening  is  purely  isotropic,  § = a and  R = a 

ij  ij 


where  o [=  J 3 J is  the  effective  stress  in  conventional  notation. 

Assuming  that  the  subsequent  yield  criterion  for  any  instan- 
taneous yield  surface  is  still  given  by,  say,  the  von  Mises  condition, 
the  yield  surface  size  parameter  k in  the  plastic  strain-rate  equa- 
tions, equations  (3.5)  and  (3.6),  is  given  by 


k = / J ?'  . §'  . 

V 2 bij  bij 


From  equations  ( f ) and  (g) 
k = — R. 

</3" 


(g) 


(h) 
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Figure  3.1  Combined  translation  and  expansion 
of  yield  surface  in  stress  space 


Figure  3.2  Schematic  of  solution  domain  and 
typical  mesh  point  types 
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In  analogy  with  the  strain  hardening  assumption  for  isotropic 
hardening,  we  assume  that  R in  equation  (h)  above  is  a function  of  the 
accumulated  plastic  strain.  In  terms  of  the  definition  of  effective 
plastic  strain 


de 


P 


This  means  that 


k = 


1 


(i) 


Differentiating  equation  (i)  with  respect  to  time,  and  noting  that 


1 p 1 p 

~ 2 £11  - ~ 2 6 ’ 


we  have  for  the  tube  problem, 


3k  _ 1 dR 
3t  ~ 3 T~P 


(j) 


Let  the  slope  of  the  universal  effective  stress/effective  plastic 
strain  curve  in  the  case  of  purely  isotropic  hardening  and  for  the 


current  accumulated  plastic  strain  be  g(eP).  That  is, 


da  _ p. 

===-  = g(e  ) • 


de1 


(k) 


Then,  for  purely  isotropic  hardening,  equation  (j)  is  rewritten  as 


3k  _ 1 
3t  _ 3 


g(eP) 


(J) 


When  there  is  combined  isotropic  and  kinematic  hardening,  let  m denote 
the  fraction  of  g(e^)  accounted  for  by  the  isotropic  part  of  the  gen- 
eral hardening  process.  We  can  then  write,  replacing  a for  isotropic 


hardening  in  equation  (k)  by  R in  the  case  of  combined  hardening, 


dR 


de 


P 


mg(e  ) . 


(m) 


Therefore,  equation  ( JL ) becomes 


Bk  1 
cit  = 3 


mg(eP) 


(3.7) 


This  is  the  equation  giving  the  timewise  variation  of  the  yield 
surface  size  parameter  in  equations  (3.5)  and  (3.6)  for  combined 
hardening. 


The  time-rate  of  change  of  the  coordinates  of  the  center  of 
the  yield  hypersurface,  in  any  dynamic  plastic  deformation  process 
described  by  purely  kinematic  hardening,  is  given  by  (see,  Goel  and 
Malvern,  1970), 


dcr  2 p 

~5t  = 3 g(e  } 


BeP 

"St" 


(n) 


BT  1 -p 

"3F  “ 3 g(e  } dF  ‘ (p) 

With  a fraction  m of  g(eP)  already  accounted  for  by  isotropic  harden- 
ing, the  remaining  fraction,  that  is,  the  part  (l-m)g(eP)  may  be 
ascribed  to  kinematic  hardening.  Thus,  we  modify  equations  (n)  and 
(p)  for  combined  hardening  as 
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P 

e in  equations  (3.7),  (3.8),  and  (3.9),  representing  current  total 
plastic  strain,  is  obtained  progressively  from  the  definition  of 
effective  plastic  strain  which,  in  the  case  of  two  non-zero  strain 
components,  reduces  to  (see  note  prior  to  equation  (j)) 


(AeP)2  + 


(3.10) 


Equations  (3.1),  (3.2),  (3.3),  (3.4),  (3.5),  (3.6),  (3.7),  (3.8), 
and  (3.9)  with  the  following  assumptions 


3e  _ 1 9ct 
3t  " E 3t  + 


1 Bt 
G 3t 


+ 


(3.11) 

(3.12) 


constitute  a complete  system  of  quasilinear  hyperbolic  partial  differ- 
ential equations,  which  can  be  solved  by  the  method  of  characteristics. 


3.4  Interior  Differential  Equations 
and  Associated  Characteristics 


We  first  substitute  equations  (3.5)  and  (3.6)  in  equations 
(3.11)  and  (3.12),  and  then  use  the  compatibility  conditions,  equations 
(3.3)  and  (3.4),  to  get 


Bu  _ 1 3a  A 
3x  _ E ciF  + E 


(3.13) 


where 


3v  _ 1 Bt  B 
^ - G cFt  + G 


E\ 


2(a-a*) 


3 


(3.14) 


exp 


-1 


(3.  15) 
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B = GX 


2 (T-T  *) 


exp  i of 


.-V  3 

C — T— 

v k 


0)- 


(3. 16) 


Equations  (3.1),  (3.2),  (3.13),  and  (3.14)  can  be  written  in  a compact 
form  as 

F w + H w = D (3. 17) 

~~t x ~ 

T A B T 

where  w and  D are  the  vectors  (u,cr,v,T)  and  (0,0,—,—)  , respectively, 

r.  It 


and  F and  H are  the  following  matrices 


F = 


p 

0 

0 

0" 

“o 

-1 

0 

0" 

0 

0 

p 

0 

0 

0 

0 

-1 

0 

1 
” E 

0 

0 

7 

H 

1 

0 

0 

0 

0 

0 

0 

- 1 

r. 

0 

0 

1 

0 

U 



The  characteristic  velocities  c for  equations  (3. 17)  are  the  roots 
of  the  equation 

Det  (c  F - H)  = 0 . 


The  determinant  is  easily  evaluated  to  obtain  the  following  four  roots 

of  c;  ± c , ± c , where 
o s 


c = 
o 


c = 
s 


'G 


There  are  four  linearly  independent  null  vectors  corresponding  to  these 
four  roots.  The  vector  for  (+CQ)  is  given  by 


l (c  F - H)  = 0 
— o — — 
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This  expands  to 


pCo  *1 


X3  = ° 


L - 


c 2 / E 
o 3 


= 0 


pCo  l2 


*4=  ° 


2„  - 


c 

o 4 


= 0 


From  the  first  two  of  these  equations  we  may  choose  2 ■=  1,  arbitrar- 
ily, giving  2 = pc  . From  the  second  two  equations,  it  is  seen  that 

o O 

2 

the  determinant  does  not  vanish,  since  pc^  4 G,  and  hence  the  only 

possible  solution  is  that  - 2^  = 0.  For  these  four  components  of 

the  null  vector  2 the  interior  differential  equation  for  c is  then 

o 

given  by 

_2T  F dw  = if  D dt 

which  simplifies  to  the  differential  relation 


dc  - p c du  = - A dt  along  dx  = + c dt  . (3.18) 

o o 

Replacing  (+cq)  by  (-c0)  for  the  negative  characteristic,  we  get 

da  + p c^  du  = - A dt  along  dx  = - cq  dt  . (3.19) 

Proceeding  similarly  for  c = +c  and  c = -c  , we  get  the  following 

s s 

two  relations,  respectively. 


dT 

- pc 

s 

dv  = 

- B dt 

along 

dx  = + c dt 
s 

(3.20) 

dT 

+ p c 
r s 

dv  = 

- B dt 

along 

dx  = - c dt  . 
s 

(3.21) 

Apart  from  the  above  four  different  relations,  there  are  five  inde- 


pendent constitutive  equations  discussed  in  Section  3.3.  For  a material 
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point,  that  is,  with  x = constant,  these  constitutive  equations  reduce 
to  ordinary,  differential  equations  (with  time  t as  the  independent 
variable)  and  thus  can  be  looked  upon  as  differential  relations  valid 
along  the  characteristic  dx  = 0.  The  constitutive  equations,  suit- 
ably simplified,  along  with  equations  (3.18)  through  (3.21),  are 
reproduced  below  for  review  as  the  final  reduction  of  the  governing 
set  of  equations  for  the  problem  under  consideration. 

Summary  of  differential  relations: 

(i)  do  - p c du  = - A dt  along  dx  = + c dt 

o o 

(ii)  do  + p c du  = - A dt  along  dx  = - c dt, 

o o 

where  A is  given  by  equation  (3.15), 


(iii) 

dT  - p c 
K s 

dv  = - B dt 

along 

dx  = + c dt 
s 

(iv) 

dT  + pc 

s 

dv  = - B dt 

along 

dx  = - c dt, 
s 

where  B is  given  by  equation  (3.16), 


(v) 

d6P 

= (A/E)  dt 

along  dx  = 0 

(vi) 

. P 

dv 

= (B/G)  dt 

along  dx  = 0 

(vii) 

dk  = 

: | g(eP)  dt \/3 (A/E) 2 + (B/G)2 

along 

dx  = 0 

(viii) 

* 

da 

= (l-m)g(eP) 

• A/E 

along 

dx  = 0 

(ix) 

* 

dT 

= | (l-m)g(eP) 

• B/G 

along 

dx  = 0 

where  eP  is  obtained  by  a progressive  integration  of  the  relation  in 
equation  (3.10)  and  g(eP)  has  the  meaning  indicated  by  equation  (k)  in 
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Section  3.3.2.  The  above  nine  equations  are  to  be  integrated,  sub- 
ject to  given  initial  and  boundary  values. 


3. 5 Initial  and  Boundary  Conditions 

Initial  and  boundary  conditions  were  chosen  from  the  1969 
report  of  Lipkin  and  Clifton.  Since  the  maximum  amount  of  experimen- 
tal results  were  reported  for  one  set  of  initial  and  boundary  data,  it 
was  decided  to  use  that  particular  set  which  is  as  follows: 


T(x,0) 
u(0 , t) 
v(0 , t) 


- 3250  psi 
= 435  in. /sec 
= 22  in. /sec 


after  a parabolic  rise  time 
of  10  )osec. 


The  boundary  values  of  velocity  steps  were  taken  with  a parabolic  rise 
only  to  prevent  any  oscillation  in  the  numerical  solutions.  The  ini- 
tial plastic  shear  strain  was  evaluated  by  first  converting  t(x,0)  to 
effective  stress,  and  finding  the  corresponding  effective  plastic 
strain  from  the  universal  stress-strain  relation  supplied  by  these 
authors  (1969) , and  then  converting  this  effective  strain  into  an 
equivalent  engineering  shear  strain.  Initial  value  of  the  yield  sur- 
face size  parameter  k,  for  the  generalized  combined  hardening  formula- 
tion, was  obtained  as  follows. 

k(x, 0)  = T (x, 0)  - T*(x,0) 

T*(x,0)  = (1-m)  [”t ( x , 0 ) - 

L S3 

where  a is  the  yield  stress  in  uniaxial  tension.  The  above  two  equa- 

y 

tions  describe  the  fact  that  the  initial  yield  configuration  is  attained 
by  a combined  expansion  and  translation  of  the  virgin  yield  due  to  a 


purely  torsional  loading. 
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3. 6 Numerical  Solution  Schemes 

Figure  3.2  (a)  shows  the  schematic  representation  of  the  domain 
in  which  computations  were  performed;  two  typical  mesh  types,  for  bound- 
ary and  interior  point  calculation,  are  also  shown  in  Figures  3.2(b) 
and  3.2  (c) . 

3.6.1  Difference  Equations  for 

Predictor-Corrector  Scheme 

Boundary  point  (Ref.  Figure  3.2  (b) ) 

* 

At  = Atc/(c+c) 
o os 

\ = (At  - A t*)/ (0.  5 At) 

XB+  UXC-V  ' 

X in  this  linear  interpolation  formula  is  a dummy  variable  and  stands 
for  N,  N , S , S , K , V,  in  turn. 

(i)  NT-Nc  + pco(UT-Uc)  = -AcAt/2 

(ii)  ST  - SR  + pcs(VT  - VR)  = - Br  At* 

(iii)  N*  - N*  = | (1-m)  gB  At(Afi/E) 

(iv)  St  " Sb  = | (1-m)  SB  At(BB/G) 

(V)  KT  - KB  = I ^B  At  / 3(Ab/E)2  + (Bg/G) 2 

(vi)  - eP  = At  J (ab/e)2  + i (bb/g)2  . 

In  the  above  equations  g stands  for  g(eB).  In  equations  (i)  - (vi)  , 

B B 

A^  implies  the  evaluating  of  the  A-function  (equation  (3.15))  with 

■if.  if 

its  five  arguments,  N,  N , S,  S , K,  representing  these  state  variables 
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at  the  point  C in  Figure  3.2  (b) , and,  similarly,  for  other  sub- 
scripted A or  B symbols.  Values  of  N^,  ST,  N , S^,  K,  e^,  obtained 
from  above,  give  the  first-order  solutions  for  these  variables. 

These  are  next  used  as  the  first  guesses  in  the  following  nonlinear 
equations  arising  out  of  a predictor-corrector  scheme  of  solution. 

We  first  obtain  the  mean  values  of  the  arguments  of  the  A and 
B functions,  that  is , N,  N , S,  S , K,  along  the  pertinent  character- 
istic, as  follows  (X  standing  for  these  variables), 

*rc  = °'5<xt  + V 
XTR  * °-5<XT  + V 

*rB-  °'5<xr  + V • 

In  the  following  Atc  will  mean  A(NTC,  Ntc>  S^,  Stc>  Ktc)  , and, 
similarly,  for  other  notations  of  this  type. 


(vii)  NT  - Nc  + pCQ(UT  - Uc)  + Atc  At/2  = 0 


(viii)  S - 


S„  + pc  (V  - V ) + B_  At  = 0 
R K s T R TR 


(ix)  N*  - N*  - | (l-m)gTB  At  (A^/E  ) = 0 
<*>  < - SB  - I 4t  <BTB/G)  = ° 


(xl)  kt  - kb  - | "«XB  4t  /=(\b/E)2  + 'BTB/G)2  = °' 


In  equations  (vii)  to  (xi) , g has  the  following  meaning: 

TB 


’TB 


P P 
®T  + 6B 


g 


2 
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Equations  (vii)  to  (xi)  constitute  a system  of  five  simultaneous 
nonlinear  algebraic  equations,  which  is  solved  by  the  Newton-Rephson 
iteration  procedure  as  discussed  in  Section  3.6.2.  If  the  iterated 
values  do  not  meet  the  desired  accuracy  conditions,  further  iteration 

p 

is  done  with  the  following  updating  of  , 

«;  = £b  + « j <\b/e>2  + 1 <btb/c>2 

with  the  latest  values  used  in  computing  A^,  , B^. 

General  interior  point  (Ref.  Figure  3.2  (c)) 

Interpolation  for  the  values  of  N,  N , S,  S , K,  V at  L and  R, 
is  done  as  follows,  X standing  for  these  variables  in  turn, 

MXA  - V 

"R  = XB  + X(XC  - V ' 

Finite  difference  forms  of  the  governing  equations  are: 

(i)  Nt  - Na  -pco(UT  - UA)  = - AA  At/2 

(ii)  Nt  - Nc  + PCq(Ut  - Uc)  = - Ac  At/2 

(iii)  ST  - SL  - pcs(VT  - VL)  = - ^ At* 

(iv)  Sx  - SR  + pcs(VT  - VR)  = - Br  At* 

(v)  N*  - N*  = | (l-m)gB  At(AB/E) 

(vi)  S*  - S*  = | (l-m)gB  At(Bfi/G) 

(vii)  kt  - Rg  = | mgB  At  J 3(Afi/E)2  + (BB/G) 2 

- 4 = At  v7  (ab/e)2  + J (bb/g)2  . 


(viii) 
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Subscripted  notations  have  the  meaning  given  in  the  discussion  of  the 
boundary  point  finite  difference  equations.  Values  of  the  state  var- 
iables obtained  from  the  above  equations  are  used  as  first  guesses  in 


solving  the  following  nonlinear  equations  originating  from  a predictor- 
corrector  format.  These  equations  are  written  after  first  eliminating 
U , VT  from  e1uations  (i) » (if),  (iff),  (iv)  given  above.  Also,  the 
arithmetic  mean  of  argument  values  between  the  top  and  base  points  in 
any  characteristic  is  taken  as  specified  before.  Thus  we  have 

= °'5<XT  + V 
*■11,  = °-5<XT  + V 

xtb=0-5<xt  + V 
XTK=  °-5<XT+  V 
*rc  = °'5<xt  + V 

where  X^,  for  example,  stands  for  the  arguments  in  the  functions  A^, 

B in  the  following: 

TB 

(ix)  ST  - 0.5[Na+Nc^co(Uc-0a)  -§  At(AIA+ATC)l  = 0 

(x)  ST  - 0.5[SL+SR+pCs(VR-VL)  -At*(BTL+BTR)]  = 0 

<xi)  N*  - NB  - | (l-m)gTB  At(ATB/E)  = 0 

<Xll)  4 - SB  - 5 (1-")gTB  4t(BTB/G)  = ° 

(xiii)  K^,  - Kr  - 1 mgTB  At  -/  3(ATB/E)  + (BTB/G)  = 0 . 

Equations  (ix)  to  (xiii)  represent  a system  of  five  nonlinear  algebraic 
equations  which  have  to  be  solved  for  N^,  S^,  N^,  S^,  K.  The  method  of 
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solution  is  discussed  next,  with  the  last  iterated  value  updating  the 
total  effective  plastic  strain  as  follows  before  a new  loop  is  computed, 

4 "b  + 4t  J <aib/e)2  + I (btb/g)2  ■ 

3.6.2  Newton-Raphson  Method  with 
Gauss-Jordan  Elimination 

Equations  (vii)  - (xi)  for  boundary  point  and  equations  (ix) 
to  (xiii)  for  an  interior  point  represent  nonlinear  equations  of  the 
form 

f1(x1'  V*’V  = 0 

(a) 

W v-'-’V  = ° 

Attempts  to  solve  these  equations  by  the  Gauss-Seidel  iterative  scheme 
did  not  prove  to  be  very  successful.  A comparison  of  the  theoretical 
solutions  so  obtained  with  the  experimental  data  indicated  that  the 
value  of  \ should  be  increased  in  the  equations  (3. 15)  and  (3. 16)  to 
have  desired  rate  insensitivity,  but  increasing  on  the  other  hand, 
required  decreasing  the  mesh  size  in  order  to  satisfy  the  sufficiency 
conditions  of  convergence  of  Gauss-Seidel  schemes  (see  Carnahan  et  al. , 
1969).  It  was  thus  decided  to  use  the  more  powerful  Newton-Raphson 
scheme.  The  solution  procedure  in  this  scheme  (Carnahan  et  al.  (1969)) 
starts  with  a guess  for  the  desired  roots,  in  this  case  obtained  from 
the  first-order  solutions  discussed  above.  Let  these  guessed  values 
be  labeled  x1Q,  x^,  x3Q , x4Q,  x , and  we  assume  that  ix  , Ax2,...,£x5 
are  the  increments  necessary  in  these  values,  respectively,  for  a set 
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of  exact  roots.  Expanding  the  functions  in  equation  (a)  about  the 
guessed  solutions  in  a Taylor  series  and  retaining  terms  up  to  the 
first  derivatives,  we  have  the  following  system  of  linear  equations, 

df . x 

(Ax  ) = - (fiQ),  i = 1 , 5 (b) 

j j = 1| . . . i 5 


where  (Bf./Bx.)  is  the  matrix  of  the  first  partial  derivatives  of  the 
i 3 

function  f.  w.r.t.  x.  (that  is,  the  determinant  of  (3f./Bx.)  repre- 
i J i 3 

sents  the  Jacobian  of  the  system  in  equation  (a)),  (Ax.)  is  the 

J 

vector  representing  increments  and  f represents  the  vector  of  the 

values  of  the  functions  f.(x  , x , ...,x  ).  The  linear  equations (b) 

l 10  JO  50 

are  solved  for  the  increments  Ax.  by  first  constructing  what  is  called 

J 

an  augmented  matrix  A as  follows. 


This  means  that  for  indices  i and  j taking  values  from  1 to  5,  matrix 
A is  of  the  order  (5  X 6) . The  Gauss-Jordan  elimination  scheme  gives 
now  an  algorithm  by  which  a matrix  B is  generated  with  one  column  less 
than  in  A,  that  is,  B is  of  the  order  (5  X 5) . The  process  actually 
reduces  the  number  of  columns  by  one  in  each  successive  step,  the 
final  transformation  of  the  original  agumented  matrix  A representing 
the  desired  roots.  Elements  of  the  matrix  B may  be  written  as  follows 
(see  James,  Smith,  Wolford,  1970), 
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b 


a.  . 
ij 


1 < i < 5 
1 < j < 6 


5,  j-1 


\A 1 < J - 6 

au  j an * ° 


(C) 


where  a and  b are  elements  of  the  old  (A)  and  new  (B)  matrix, 

respectively.  Repeated  application  of  equations  (c)  till  matrix  B 

is  of  the  order  (5  x 1)  gives  the  roots  AX.,  j = 1,...,5. 

J 

In  the  spirit  of  the  Newton-Raphson  procedure,  the  refinement 
on  the  guessed  roots  is  given  by 


x.  . = x.  + Ax.  , i = 1 , . . . ,5  . (d) 

1,1  1,0  1 

If  these  newly  obtained  roots  do  not  satisfy  equations  (a),  then  the 
process  is  repeated.  The  following  tolerance  limit  is  reached  at  exit 


. 999 (N) < 


I 


i=l 


X. 

l ,r+l 
X. 

i,r 


<1.001  (N) 


(e) 


where  N is  the  number  of  variables  involved  and  (r+1)  is  the  number 
of  iterations. 

In  the  computations,  the  numerical  solution  of  equations  (a) 
was  hindered  by  a systematic  oscillation  between  two  successive  iter- 
ations, apparently  caused  by  an  inflection  point  near  the  desired  root 
This  problem  was  overcome  by  a system  of  fractional  incrementing  such 
that  equation  (d)  is  replaced  by 


where  3 


x.  = 

i ,r+l 

is  a fraction. 


x.  +3  Ax.  (f) 

i , r i 

This  scheme  was  made  operative  only  after 


five  iterations,  if  the  desired  accuracy  was  not  obtained  till  then. 
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In  any  case,  computations  at  a given  mesh  point  were  allowed  up  to 
a maximum  number  of  15  iterations.  In  all  cases  where  this  occurred, 
the  final  error  was  greatly  reduced  compared  to  that  with  first- 
order  solutions. 

3.7  Results  - Comparison  with  Experimental 
Data  and  Rate-Independent  Solution 

Results  of  the  numerical  solutions  are  discussed  in  the  follow- 
ing, in  the  light  of  the  experimental  data  and  a rate-independent 
analysis  given  by  Lipkin  and  Clifton  (1969),  also  by  Lipkin  (1971). 

Figures  3.3  through  3.10  present  the  results  for  the  isotropic 
hardening  assumption.  The  numerical  solution  represented  by  them  was 
obtained  with  a step  size  of  AT  = 0.5  psec  for  the  values  of  \ = 300 
and  a = 10  in  equations  (3.15)  and  (3.16).  The  parameters  X and  a 
control  the  degree  of  rate-sensitivity,  with  higher  values  for  them 
indicating  less  rate-dependence,  in  these  equations.  The  extent  of 
rate-insensitivity  implied  by  the  numbers  used  here  for  X and  a can  be 
realized  by  comparing  the  rate-dependent  solution  with  the  rate- 
independent  one. 

Figure  3.3  shows  the  e - t and  Ay  - t plots  for  X = 5.75  inches. 
For  the  longitudinal  strain,  the  constant  state  region  of  rate-inde- 
pendent theory  following  the  fast  wave  region  matches  well  with  exper- 
iments up  to  about  87  psec,  and  this  result  is  reproduced  by  the  rate- 
dependent  solution  as  well.  The  experiments  show  a perceptible  rate 
of  straining  beyond  87  psec,  while  the  rate-independent  solution  con- 
tinues to  predict  a constant  state  up  to  120  psec.  As  suggested  by 
Lipkin  and  Clifton  (1969) , a rate-dependent  theory  such  as  the  present 


Rate-dependent  solution 
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Figure  3.3  (-e)-t  and  Ay-t  plots  for  x = 5.75  inches 
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Figure  3.4  Comparison  of  theoretical  and  experimental  strain  trajectories 
at  x = 5.75  inches 
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one  is  inherently  capable  to  explaining  this  non-constant  state,  though 
in  the  present  situation  the  experimental  response  is  actually  lagging 
the  rate-dependent  solution.  With  the  shear-strain  record,  the  rate- 
dependent  solution  almost  reproduces  the  fast  wave  region  solution  of 
rate-independent  theory,  but  again  the  experimental  results  lag  the 
theoretical  solutions.  The  rate-dependent  solution  stops  showing  the 
approximate  constant  state  almost  at  the  same  instant  with  the  exper- 
iment, while  the  rate-independent  solution  continues  to  yield  a con- 
stant state  up  to  about  120  psec.  In  general,  the  experimental  response 
again  lags  both  the  theoretical  solutions  in  the  non-constant  state 
region.  Finally,  it  is  seen  from  the  trend  of  the  curves  at  the  instant 
up  to  which  computations  were  done,  that  the  constant  state  region 
following  the  slow  simple  waves  of  rate-independent  theory  came  much 
too  soon;  also  the  maximum  values  of  e and  Ay  obtained  by  the  rate- 
dependent  solution  are  closer  to  the  experimental  maxima  than  the  rate- 
independent  solution. 

Figure  3.4  shows  the  strain  trajectories  corresponding  to 
Figure  3.3.  The  rate-dependent  solution  presented  is  seen  in  general 
to  be  closer  to  the  experimental  response  than  the  rate-independent 
solution. 

Figure  3.5  shows  the  e-t  plots  at  four  stations,  x = 0, 
x = 1.875  inch,  x = 3.375  inches  and  x = 5.75  inches.  (Experimental 
data  and  rate-independent  solution  for  the  first  three  of  these  sta- 
tions are  not  available. ) The  rate-dependent  solutions  are  seen  to 
reproduce  the  feature  of  a rate-independent  solution  in  that  the  con- 
stant state  duration  increases  with  increasing  distance  from  the  bound- 
ary. Similar  observations  can  be  made  from  Ay  - t plots  presented  in 


75 


Figure  3.6  for  these  four  stations.  In  both  of  these  figures  it  is 
seen  that  the  constant  velocity  inputs  translate  themselves  into  con- 
stant strain  outputs  at  the  boundary  with  different  rise  times  for  e 
and  Ay  from  that  for  the  velocity.  Also,  the  interior  strains  are 
seen  to  take  progressively  longer  time  to  reach  their  respective  plateaus. 

Figures 3. 7 (a)  and  (b)  show  the  velocities  of  propagation  for 
different  levels  of  longitudinal  strain.  Figure  3.7  (a)  Considers  the 
strain  levels  corresponding  to  a fast  simple  wave  solution  in  rate- 
independent  theory.  It  is  seen  that  the  rate-dependent  theory  repro- 
duces better  the  propagation  speeds  between  stations  2 and  3,  i.e. , 
between  x = 3.375  inches  and  x = 5.75  inches,  while  in  general  the 
rate-dependent  solution  predicts  faster  speeds  than  the  rate-independ- 
ent theory.  More  significantly,  it  is  seen  that  the  present  rate- 
dependent  solution  shows  an  approximately  constant  speed  of  propaga- 
tion for  any  given  level  of  strain,  much  like  a rate-independent  theory. 
The  experiments  however  do  not  show  this  constancy. 

In  Figure  3.7  (b) , where  the  strains  for  the  slow  wave  region 
of  rate-independent  theory  are  presented,  the  situation  is  made  clearer 
with  the  rate-independent  solution  forming  the  lower  bound  and  the 
rate-dependent  solution  an  upper  bound  for  the  experimental  results. 

The  results  for  stations  2 and  3 are  better  represented  by  a rate- 
dependent  theory,  while  data  for  stations  1 and  2 are  closer  to  the 
rate-independent  solution. 

Again,  the  rate-dependent  solution  gives  approximately  a con- 
stant speed  of  propagation  for  any  given  level  of  strain,  unlike 
the  experiments.  In  general,  these  curves  show  that  the  rate- 
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Figure  3.5  (-e)-t  plots  for  x = 0,  1.875,  3.375,  5.75  inches 
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Figure  3.6  Av-t  plots  for  x = 0,  1.875,  3.375,  5.75  inches 


78 


i 


to 

■o 

o 


in 

o 

o 


o 

o 


o 

0 

CA 

\ 

1 *■ 

X 

o 

cn1 

« 

in 

I 

1 

CO 

c 

• 

• 

• 

• 

CO 

p 

O 

3 

- 3 

0 

-2 

o 

o 

o 

l-° 

o 

(O 

o 

■o 

o 


fl 

o 


p 

o 

a 

o 

-P 

0 


co 

i 


P 

1 — ( 

Cl) 

| 

P 

o 

03 

P 

tfl 

1 

o 

P 

COO 

to 

-p 

Cl 

-p 

s 

0 

T3 

§ 

0) 

Cl 

■a 

0 

p 

a 

ft 

0 

o 

0 

CO 

a 

-a 

I 

X 

o 

i 

03 

>i 

03 

O 

N-X 

■a 

0 

P 

c 

-P 

ca 

Jh 

•H 

ca 

£ 

03 

w 

i 

p 

0 

P 

0 

•H 

P 

•p 

e 

P 

£ 

o 

03 

•H 

p 

p 

P 

CO 

<H 

CA 

s 

P 

Q 

. . 

0 

p 

co 

o 

T 3 

<H 

i 

C 

CO 

* 

0 

II 

bJD 

«. 

0 

l— l 

Q 

« 

in 

t>  r 

oo  m 

CO  r4 


co  co 

i i 
i— i i— t 


<•) 

I 


o 

,o 

o 


w 

a 

o 

P 

P 

a 

P 


CA 

£ 

o 

•H 

4-> 

08 

P 


m 

co 

03* 


CO 

I 

03 

CA 

P 

O 

•H 

P 

03 


03  CO  CO  CO 


P 

G 

0) 

£ 


p 

£ 

0) 

£ 


P 

C 

0 

S 


lO 

I 


o 

0 

ca 


X 


X 

o 

c 

•H 


p 

p 

p 

0 

0 

.0 

ft 

ft 

ft 

X 

X 

X 

w 

w 

w 

o 

X 

• 

Figure  3.7  Velocity  of  propagation  for  various  levels  of  longitudinal  strain 
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Figure  3.8  Velocity  of  propagation  for  various  levels  of  incremental  shear  strain 
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dependent  solutions  actually  would  have  fared  better  with  still  more 
rate  insensitivity  incorporated  in  the  governing  equations,  i.  e.  , using 
still  higher  values  of  \ and  o'  in  equations  (3.15)  and  (3.16). 

Figures  3.8  (a)  and  (b)  present  results  corresponding  to 
Figures  3.7  (a)  and  (b)  for  incremental  shear  strain.  For  the  strains 
in  the  fast  wave  region  of  rate-independent  theory,  both  the  theories 
predict  higher  propagation  speeds  than  in  the  experiment.  Strains 
corresponding  to  the  slow  wave  region  of  rate-independent  theory  are 
in  general  better  represented  by  the  rate-dependent  solution  given. 

Figure  3.9  shows  the  plots  of  and  k against  time  for  the 

four  stations  specified  before.  The  difference  in  the  ordinates  for 
any  pair  of  curves,  i.e.  , a/j^  - k for  any  station,  represents  the  over- 
stress and  the  smallness  of  this  quantity  shows  the  large  extent  of 
the  rate-insensitivity  incorporated  in  the  present  rate-dependent  solu- 
tion. These  solutions  also  show  a continual  increase  of  stress  at 
the  interior  stations  signifying  "total  loading”  (Cristescu,  1967), 
while  at  the  boundary  there  was  even  a dynamic  relaxation. 

Figure  3.10  shows  for  the  four  stations  the  stress  paths  in 
the  ct-T  plane,  during  the  dynamic  plastic  deformation  caused  by  the 
stress  waves  passing  through  them.  What  is  remarkable  here  is  that, 
excepting  at  the  boundary,  the  stress  trajectories  look  very  much  like 
those  for  a rate-independent  solution  for  step-loading  of  a pre- 
twisted tube  (Clifton,  1966).  In  another  sense,  this  shows  the  extent 
of  rate-insensitivity  represented  by  the  rate-dependent  solution  given. 

Figure  3.11  presents  the  strain  time  plots  for  x = 5.75  inches 


under  the  assumptions  of  kinematic  hardening,  while  Figure  3.12  presents 
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the  corresponding  results  for  combined  kinematic  and  isotropic  harden- 
ing. Both  of  these  results  were  obtained  for  the  step  size  AT = 1 psec 

-4  -5 

with  X = 10  and  a = 6 for  kinematic  hardening,  and  X = 10  and 

O' =10  for  combined  hardening.  Higher  values  for  X and  ot  than  these 
in  each  case  were  causing  extremely  slow  or  oscillatory  convergence  of 
the  solutions  at  various  mesh  points.  The  agreement  between  either  of 
these  solutions  with  the  experiments  is  apparently  very  poor.  In  both 
Figures  3.11  and  3.12,  are  plotted  the  solution  for  isotropic  harden- 
ing assumption  when  the  step  size  was  AT = 1 psec.  This  solution  incor- 
porates almost  the  threshold  values  of  X and  a for  this  step  size 
(X=l  and  a=  10),  beyond  which  convergence  was  not  obtained  even  with 
15  Newton-Raphson  iterations.  The  subsequent  improvement  in  the  results 
for  isotropic  hardening  after  halving  the  step  size  (that  permitted  the 
use  of  values  of  \=  300  with  <y=  10)  shows  that  the  results  for  kine- 
matic and  combined  hardening  could  also  be  improved.  However,  the 

computer  time,  even  for  isotropic  hardening  that  involves  two  equa- 

* * 

tions  less  to  solve  at  each  time  with  o = T = 0,  is  about  13  min- 
utes on  an  IBM  360/65  machine  for  the  reduced  step  size  of  0.5  psec. 

The  kinematic  and  combined  hardening  solutions  were  therefore  not 
pursued  any  further,  also  for  the  additional  reason  that,  for  the  same 
step  size,  the  isotropic  hardening  solution  was  comparatively  better. 

Figures  3.13  and  3.14  finally  show  the  stress  trajectories, 
corresponding  to  the  strain-time  plots  of  Figures  3.11  and  3.12,  for 
the  kinematic  and  combined  hardening  assumptions,  respectively.  The 
paths  followed  by  the  center  of  the  yield  hypersurface  are  also  noted 
in  these  figures.  These  figures  show  the  large  extent  of  rate- 
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Figure  3.11  Comparison  of  strain-time  plots  for  x = 5.75  inches  with 
kinematic  hardening  and  isotropic  hardening  assumptions 
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Figure  3.12  Comparison  of  strain-time  plots  for  x = 5.75  inches  with 
combined  hardening  and  isotropic  hardening  assumptions 


Legend: 

A = Trajectory  of  yield  surface  center 
B = Initial  yield  surface 


o C = Final  yield  surface 

(psi) 


Figure  3.13  Stress  trajectory  at  x - 5.75  inches 
under  kinematic  hardening  assumption 
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Legend: 

A = Trajectory  of  yield  surface  center 
B = Initial  yield  surface 
C = Final  yield  surface 


Figure  3.14  Stress  trajectory  at  x = 5.75  inches 
under  combined  hardening  assumption 
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dependence  implied  by  the  solutions  presented.  Considering  that  the 
physical  problem  involved  an  essentially  rate-insensitive  material, 
the  disagreement  between  the  experiment  and  the  solutions  for  kine- 
matic and  combined  hardening  is  thus  only  to  be  expected  for  the 
particular  choices  of  X and  o'  in  these  solutions. 

3.8  Conclusions 

1.  It  is  shown  that  for  a material  whose  dynamic  response  can 
be  fairly  well  represented  by  a rate-independent  theory,  except  for  a 
few  details,  a rate-dependent  theory  can  also  be  made  to  give  a reason- 
ably good  agreement  with  a suitable  choice  of  its  rate-sensitivity 
parameters.  In  the  particular  case  of  the  experimental  results  anal- 
yzed here,  the  rate-dependent  theory  actually  explains  better  the  non- 
constant state  regions  in  the  experiment,  while  retaining  the  other 
good  features  of  the  rate-independent  solution. 

2.  The  difficulty  of  representing  an  essentially  rate- 
insensitive  material  response  by  means  of  a rate-dependent  theory  is 
one  of  poor  or  oscillatory  convergence  of  the  equations,  even  with  the 
Newton-Raphson  schemes  that  are  known  to  have  superior  convergence 
properties  than  the  conventional  iterations  like  the  Gauss-Seidel 
method. 

3.  In  view  of  the  fact  that  fundamental  quantities  like  the 
speed  of  propagation  of  any  given  level  of  strain  are  not  well  repre- 
sented by  either  a rate-independent  theory  or  a rate-dependent  theory, 
we  emphasize  the  need  for  further  experimental  verification,  with 
strain  measurements  done  preferably  by  diffraction  gratings  (Bell , 

1968)  rather  than  strain  gauges  which  can  show  a delayed  response 


to  large  amplitude  strain. 


CHAPTER  4 


UNLOADING  CONSIDERATIONS  FOR  COMBINED  STRESS  WAVES 

In  this  chapter  we  apply  the  rate-dependent  theory  formulated 
in  the  previous  chapter  to  some  unloading  problems  in  combined  stress 
wave  propagation.  An  important  problem  in  this  is  the  determination 
of  the  loading/unloading  boundary  also  called  the  elastic/plastic 
boundary,  and  before  we  discuss  the  solutions  of  some  specific  prob- 
lems, certain  essential  features  connected  with  this  boundary  will  be 
discussed  in  the  next  two  sections. 

4.1  Equations  Governing  Elastic/Plastic  Domains 

By  definition,  the  loading/unloading  boundary  is  a moving 
boundary  in  the  physical  body  represented  by  a fixed  line  in  the 
characteristic  plane  that  separates  regions  in  which  the  material 
response  is  elastic  from  those  in  which  it  is  plastic.  These  moving 
elastic/plastic  boundaries  are  called  unloading  (loading)  waves  if 
the  material  at  a section  changes  from  a plastic  (elastic)  state  to 
an  elastic  (plastic)  state  as  the  wave  passes  the  section. 

It  is  thus  clear  that  the  equations  in  the  plastic  or  loading 
domain  remain  the  same  as  given  in  Chapter  3.  The  equations  govern- 
ing material  response  in  the  elastic  domains  are  obtained  by  simply 
setting  the  functions  giving  the  plastic  strain  rate,  i.e. , functions 
A and  B of  equations  (3.15)  and  (3.16)  to  zero  in  all  the  governing 
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equations.  The  consequence  of  this  is  that  a , T , k,  e , e , and  y 
remain  constant,  at  the  respective  maximum  values,  for  any  particular 
x undergoing  unloading. 

The  criterion  to  decide  whether  to  use  the  plastic  or  elastic 

equations  at  a particular  mesh  point  depends,  for  the  rate-dependent 

< 2 

theory  followed  in  this  study,  on  the  inequality  J > k where  the 

stresses  in  the  expression  for  J are  computed  as  though  the  govern- 

ing  equations  are  elastic  and  k is  the  maximum  value  of  work  hardening 

2 

for  the  section.  If  J < k , the  equations  governing  are  elastic,  and 
the  domain  where  this  holds  is  the  elastic  or  unloading  domain. 

4. 2 The  Initial  Slope  Problem 

The  determination  of  the  slope  of  the  loading/unloading  bound- 
ary near  the  impacted  end,  at  the  instant  when  the  applied  load  begins 
to  decrease,  can  be  a difficult  problem  in  rate-independent  theory,  as 
shown  by  the  works  of  Clifton  (1968)  and  Ting  (1969,1972).  In  the 
following  we  use  the  stress  conditions  of  Ting  (1969)  that  would  give 
unloading  of  a specific  character  in  a rate-independent  material  and 
compare  their  effects  on  materials  exhibiting  a strain-rate  effect. 

4.2.1  Comparison  with  the  Rate- 
Independent  Situation 

Ting  (1969) , in  considering  several  cases  of  unloading  from 
the  impacted  end  of  a semi-infinite  tube,  discusses  the  situation 
where  the  material  is  already  in  a plastic  state  at  the  boundary  when 
unloading  occurs.  Two  subclasses  of  relevance  to  the  present  study  are 
described  by  the  following  conditions  on  the  stress  derivatives  before 
and  after  the  instant  of  unloading,  at  the  boundary  x = 0: 
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Case  (i) 


6 

c 


Case  (ii) 


a3.  < 6 

t c 


with 


where 


Subscript  t indicates  the  time  derivative,  superscripts  b and  a refer 

to  these  derivatives  "before”  and  "after"  the  instant  of  unloading, 

and  c is  the  "slow"  wave  speed  of  rate-independent  theory  (see 
s 

Clifton,  1966). 

In  the  present  study  a numerical  investigation  was  made  of 
the  above  two  cases.  Ting  (1969)  bases  his  rate-independent  theory 
on  the  assumption  of  isotropic  hardening  and,  for  the  sake  of  compar- 
ison, the  same  assumption  is  made  in  the  rate-dependent  formulation. 
The  quantity  c^  in  the  above  inequalities  was  calculated  from  the 
formula  given  by  Clifton  (1966).  The  stress  pulses  in  ct  and  T at 
the  boundary  were  so  constructed  as  to  satisfy  these  inequalities  for 
the  two  cases  considered.  This  was  achieved  by  a linear  rise  of  a 
and  T (at  different  rates)  and  a linear  (in  some  cases,  almost  sudden) 


C(0 , t)  CT(0,t) 

T (0  , t)  T (0  , t) 
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Figure  4.1  Unloading  for  Case  (i)  with 
no  dwell  period 
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o co 
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Figure  4.2  Unloading  for  Case  (ii)  with 
no  dwell  period 
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decay  in  c combined  with  a linear  decay  in  T.  Two  pulse  shapes  were 
considered  for  each  of  Cases  (i)  and  (ii)  , one  with  a dwell  period, 


(OOt-°)>  and  one  without,  for  the  maximum  values 

attained  in  a and  T.  The  computer  program  used  was  for  loading  and 
unloading  in  a semiinfinite  tube,  with  a time  step  of  1 psec.  The 
logic  for  finding  whether  a certain  mesh  point  is  within  the  loading 
or  unloading  domain  is  as  stated  in  Section  4.1.  Rate  sensitivity 
parameters  used  were  X = 1,  a = 10  in  equations  (3.15)  and  (3.16). 


. 4.2.2  Corresponding  Results  and  Conclusion 

Figures  4.1  and  4.2  represent  the  results  for  Cases  (i) 

( \b  b 

and  (ii)  , respectively,  with  > 0 , cr^_  0.  Details  about  the 

pulse  shapes,  conforming  with  the  conditions  for  these  two  cases,  are 

also  shown  in  these  figures.  In  the  characteristic  plane,  the  solid 

line  in  each  figure  is  obtained  by  joining  the  outermost  unloading 

points  for  the  first  few  rows.  It  is  seen  that  the  unloading  wave 

speed  c in  Case  (i) , Figure  4.1,  is  less  than  c , the  elastic  shear 
u.  ^ 

wave  speed.  For  this  same  case,  Ting's  rate-independent  theory  pre- 
dicts c > c„.  Of  additional  contrast  is  the  fact  that  unloading  does 
u 2 

not  start  at  the  same  instant  that  the  stresses  begin  to  decrease  at 
the  boundary.  The  computer  solution  based  on  the  rate-dependent  theory 
had  shown  that  there  was  a dynamic  relaxation  involved  at  the  boundary, 
past  the  initial  instant  of  decrease  for  the  applied  stresses,  and 
elastic  unloading  only  began  when  the  relaxation  boundary  was  reached. 

Figure  4.2  shows  the  results  for  Gtese  (ii)  , and  it  is  seen  that 
c^  < c^,  in  agreement  with  Ting's  conclusion  for  this  case.  Also,  since 
the  solid  line  joins  the  outermost  unloading  points  in  different  rows, 


t (^sec) 


t 


Figure  4.3  Unloading  for  Case  (i) 
with  dwell  period 


t 


Figure  4.4  Unloading  for  Case  ( i i ) 
with  dwell  period 
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unloading  may  be  said  to  have  occurred  in  this  case  at  almost  the 
same  instant  when  stresses  began  to  decrease  at  the  boundary.  For  the 
particular  pulse  geometry  considered,  the  computer  solution  showed 
a rapid  dynamic  relaxation,  and  this  is  reflected  in  the  almost  instan- 
taneous unloading  response  of  the  material. 

Figures  4.3  and  4.4  show  the  results  of  the  same  two  cases, 
but  with  a dwell  period  in  the  stress  pulses.  Figure  4.3  shows  that, 
unlike  in  Figure  4.1,  cu  > c2’  as  Tins' s theory  predicts  for  this  case, 
though  there  is  a slightly  delayed  response  to  the  imposed  decrease  of 
the  boundary  stresses.  It  may  be  concluded  that  the  dwell  period  had 
the  effect  of  bringing  the  stresses  essentially  close  to  the  relaxation 
boundary  and  thus  material  response  in  unloading  was  similar  to  that 
given  by  a rate-independent  theory. 

Figure  4.4  shows  the  effects  of  a dwell  period  in  the  condi- 
tions of  Case  (ii).  The  effect  of  unloading  is  so  severe  in  this  case 
that  the  result  is  an  elastic  shock  wave  travelling  at  c . The  only 
difference  between  the  two  Case  (ii)  problems  represented  in  Figures 
4.2  and  4.4  is  that  in  Figure  4.4  the  dwell  period  had  served  to  make 
k,  the  yield  size  parameter,  bigger  and  this  made  the  imposed  decrease 
at  the  boundary  in  the  value  of  v/J~  more  pronouncedly  below  k. 

/ft 

In  summing  up,  one  may  generally  say  that  the  unloading  response 
of  a rate-dependent  material  is  quite  different  from  that  of  a rate- 
independent  material.  Because  of  the  dynamic  relaxation  implied  by 
rate-dependent  theory  unloading  response  may  be  delayed  with  respect 
to  the  actual  decrease  of  the  boundary  stresses.  If,  however,  the 
loading  is  of  such  a duration  that  the  stresses  reach  the  relaxation 
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boundary  before  the  decrease  in  boundary  stresses  occurs,  then 
response  of  a rate-sensitive  material  is  essentially  similar  to  that 
of  a rate-insensitive  material,  provided  of  course  the  imposed 
decrease  at  the  boundary  is  not  too  sudden. 

4. 3 Unloading  from  the  Impact  End 

of  a Prestressed  Semiinfinite  Tube 

In  this  section  we  consider  the  combined  stress  unloading 
problem  in  a pretorqued  semiinfinite  tube.  Boundary  values  prescribed 
are  of  the  mixed  type,  velocities  specified  during  loading  and  stresses 
during  unloading.  It  may  be  mentioned  that  such  information  can  be 
conveniently  extracted  from  typical  experimental  arrangements,  and  has 
actually  been  used  recently  by  Hsu  (1972)  from  his  experiments. 


4.3.1  Different  Pulse  Shapes  Considered 
and  Numerical  Scheme 

The  boundary  loading  was  represented  by  a longitudinal 
velocity  of  input,  while  the  tangential  velocity  v was  held  at  zero. 
In  three  cases  computed,  the  velocity  U was  maintained  constant  for 
a dwell  period  of  120  psec,  60  psec,  and  zero  (i.e.  , no  dwell),  in 
each  case  after  a parabolic  rise  time  of  10  psec.  Unloading  in  each 
case  was  represented  by  the  compressive  stress  a,  showing  a parabolic 
decay  over  a period  of  10  psec.  The  maximum  value  of  a from  which 
the  decrease  starts  is  determined  from  the  conditions  operative  at 
the  end  of  the  loading  era.  Results  for  the  case  with  60  psec  dwell 
period  during  loading  were  qualitatively  similar  to  those  for  the  one 
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Figure  4.5  Mixed  boundary  conditions  considered 
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with  120  psec  dwell  period,  and  are  not  presented  here.  Figures 
4.5  (a)  and  {b)  explain  the  boundary  conditions  used  in  the  computa- 
tions. 

The  numerical  scheme  followed  is  essentially  the  same  as  that 
for  Chapter  3 in  the  case  of  isotropic  hardening.  One  modification 
was  made,  however,  for  conveniently  recognizing  a point  in  the  elastic 
domain  of  the  loading/unloading  boundary.  This  consisted  in  doing 
the  first  iteration  based  on  the  elastic  equations,  in  all  cases,  and 
checking  whether  computed  from  stresses  so  calculated  is  less  than 

k at  the  previous  time  step  or  not.  If  the  answer  is  yes,  then  compu- 
tations for  the  mesh  point  in  question  are  accepted  as  final;  if  the 
answer  is  no,  then  the  plasticity  equations  are  used  as  usual.  This 
technique  for  computing  an  unloading  point  is  due  to  Cristescu  (1971). 
In  the  absence  of  any  unloading  experimental  data  to  be  matched,  the 
present  theoretical  solution  was  accepted  with  more  rate-dependence 
manifested  by  it  than  the  solution  presented  in  Chapter  3.  The  rate- 
sensitivity  parameteis  used  were  \ = 1,  a - 10.  Step  size  of 
AT  = 1 psec  was  good  for  convergence. 

4.3.2  Results 

Figures  4.6,  4.7,  and  4.8  present  the  results  of  the 
numerical  computation  for  the  boundary  conditions  given  in  Figure 
4.5  (a)  with  a dwell  period  of  120  psec.  Prestress  in  the  specimen 
is  taken  to  be  3250  psi  in  shear. 

Figure  4.6  gives  plots  of  V/J9  and  k at  x = 0,  1.88  inches 
and  5.75  inches.  Before  the  stresses  decrease  at  the  boundary  at 
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Figure  4.6  Values  of  and  k at  x = 0,  1.88,  5.75  inches  for  the 
pulse  of  Figure  4.5  (a)  acting  on  a semiinfinite  tube 
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130  psec,  the  difference  between  and  k is  seen  to  be  significant, 

and  this  represents  a fairly  rate-sensitive  material.  It  is  seen  that, 
due  to  this  rate  effect,  elastic  unloading  at  the  boundary  starts  about 
4 psec  after  the  stresses  start  decreasing  at  the  boundary.  Also,  from 
the  plots,  it  may  be  inf  erred,  that  the  equivalent  stress  V3.T  starts 
increasing  slowly,  once  unloading  takes  place  through  a rapid  parabolic 
decay. 

Figure  4.7  presents  the  strain-time  plots  at  x = 1.88  inches 
and  x = 5.75  inches.  The  effect  of  the  unloading  can  be  qualitatively 
judged  from  the  curves  for  x = 5.75  inches  and  comparing  them  with  the 
isotropic  hardening  solution  presented  in  Figure  3.11.  These  two  sets 
of  curves  were  obtained  with  identical  prestress,  rate-sensitivity 
parameters  and  step  size,  and  the  only  differences  were  a small  tan- 
gential velocity  at  the  boundary,  V = 22  ft/sec  for  the  case  repre- 
sented by  Figure  3.11,  and  the  imposed  boundary  stress  decay  for  the 
case  in  Figure  4.7.  It  is  seen  that  the  longitudinal  strains  at  the 
interior  almost  reach  their  plateaus,  while  the  shear  strains  continue 
to  rise  after  unloading  occurs.  This  situation  is  understood  by 
looking  at  Figure  4. 8 where  plots  are  given  of  a and  T versus  t for 
x = 0,  1.88  inches  and  5.75  inches.  It  is  seen  that,  as  a consequence 
of  unloading,  there  is  a larger  decrease  in  (-c),  and  hence  in  (-ct/E), 
followed  by  its  negligible^ increase,  while  T and  hence  AT/G  continues 
to  increase  and  actually  speeds  up.  Thus  the  Ay  - t curve  is  essenti- 
ally not  interrupted  by  the  unloading  wave. 
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Figure  4.7  Strain-time  plots  at  x = 0,  1.88,  5.75  inches  ior 
unloading  in  a semiinfinite  tube 
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Figure  4.8  (-cr)-t  and  T-t  plots  at  x = 0,  1.88,  5.75  inches  for 

unloading  in  a semiinf inite  tube 
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Figure  4.9  Results  of  unloading  in  a semiinfinite  tube 
for  the  pulse  of  Figure  4.5  (b) 
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Figure  4.9  presents  the  results  for  the  same  prestress,  but 
with  the  boundary  conditions  given  in  Figure  4.5  (b) . In  Figure  4.9  (a) 
are  shown  the  boundary  behavior  of  a,  T,  ,Jj~  and  k against  time.  The 
resultant  strain-time  response  at  x = 1.88  inches  and  x = 5.75  inches 
are  given  in  Figure  4.9  (b) . The  continued  rise  in  shear  strains  after 
unloading  can  be  explained  similarly  as  above. 

4.4  Unloading  in  a Pretwisted  Finite  Tube 
• Due  to  Load  Decay  in  Impact  End,  with 
Other  End  Free  to  Axial  Motion 

The  combined  effects  of  unloading  in  a finite  tube  due  to  load 
decay  on  one  end  and  a free  end  at  the  other  is  studied  next.  The 
"free"  end  in  this  case  is  actually  constrained  against  rotation  (in 
order  to  make  the  pretwist  possible) , but  is  free  to  axial  motion 
caused  by  an  impact. 

4.4.1  Pulse  Shapes  and  Prestress 
Levels  Considered 

To  see  the  effects  of  the  unloading  from  the  impact  end  and 
the  free  end,  first  separately  and  then  together,  it  was  decided  to 
choose  a dwell  period  for  the  loading  pulse  such  that  unloading  occurs 
at  the  impact  boundary  first  due  to  the  load  decay  itself,  before  the 
first  unloading  wave  reaches  it  from  the  free  end.  The  tube  length 
was  taken  such  that  the  first  unloading  wave  returns  at  100  psec,  and 

the  dwell  period  in  Figure  4.5  (a)  was  taken  as  60  psec.  In  another 

% 

problem  solved  in  this  category,  the  boundary  condition  of  Figure  4.5 
(b)  was  maintained.  Regarding  the  prestress  level,  computations  were 
done  with  two  values  of  shear  stress,  3250  psi  and  6000  psi.  It  was 
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seen  that  this  changes  the  stress-strain  distribution  in  the  tube 
both  qualitatively  and  quantitatively.  Results  are  presented  for 
6000  psi  prestress  for  a loading  pulse  of  60  psec  dwell  period  in 
Figure  4.5  (a),  and  3250  psi  for  the  boundary  condition  of  Figure 
4.5  (b). 

4.4.2  Loading/Unloading  Domains 
and  Other  Results 

Figures  4.10  through  4.14  refer  to  the  problem  with  6000  psi 
prestress  in  shear  and  60  psec  dwell  period  in  the  boundary  condition 
of  Figure  4. 5 (a) . 

Figure  4.10  shows  on  the  characteristic  plane  the  isolated 
domains  of  loading  and  unloading  obtained  for  this  problem.  The 
demarcation  boundaries  separating  loading  and  unloading  regions  were 
actually  obtained  by  joining  the  points  satisfying  the  unloading 
criterion.  Locations  of  two  interior  stations  (marked  I and  II)  at 
which  stresses  and  strains  are  reported  in  the  following  are  also 
shown.  Actually,  the  characteristic  diagram  shows,  for  example,  when 
such  stations  were  undergoing  elastic/plastic  loading  and  when  they 
were  in  the  process  of  elastic  unloading. 

Figure  4.11  presents  plots  of  */j7  and  k for  the  two  boundaries 
and  the  two  interior  stations.  It  is  seen  that  the  location  of  the 
second  station  is  such  that  it  undergoes  unloading  due  to  reflection 
from  the  free  end  even  before  the  impact  boundary  experiences  unload- 
ing. The  free  end  is  of  course  always  an  unloading  point.  Up  to  the 
maximum  time  for  the  impact  boundary  for  which  computations  are  reported 
here,  station  II  exhibits  plastic  reloading  twice,  while  station  I has 


one  plastic  reloading  and  one  elastic  reloading. 
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Figure  4.11  a/J,  k plots  for  the  finite  tube  at 
9.915  inches  (right  boundary) 
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Detailed  behavior  of  the  stresses  a and  T at  the  two  boundaries 

and  the  two  interior  stations,  not  easily  visualized  by  intuition,  are 
given  in  Figure  4.12.  It  is  seen  that  the  longitudinal  stresses  at 

stations  I and  II  changed  their  nature  from  compressive  to  tensile,  which 

in  a repeated  cycling  means  that  a vibratory  process  was  going  on  dur- 
ing which  the  permanent  strain  was  increasing  and  decreasing  period- 
ically. Shear  stresses  at  the  two  boundaries  are  also  seen  to  change 
their  directions  at  least  once  before  the  final  time  of  computation. 

Figure  4.13  shows  the  e-t  and  Ay  - t plots  for  the  two  bound- 
aries and  the  two  inner  stations.  Of  these  the  e-curve  for  the  impact 
boundary  and  the  Ay-curve  for  the  free  end  are  directly  understood 

from  the  a and  T curves  for  the  impact  and  free  end,  respectively,  in 

Figure  4.12.  The  other  curves  are  not  so  easily  understood  and  must 
be  taken  as  the  result  of  multiple  interactions. 

Figure  4.14  finally  presents  the  stress-trajectories  at  the 
two  interior  stations.  Information  like  the  change  of  stress  from 
a compressive  to  tensile  character  and  the  number  of  times  plastic 
reloading  occurs  are  clearly  brought  out  in  this  figure. 

Figure  4.15  shows  the  results  for  a prestress  in  shear  of 
3250  psi  and  the  boundary  condition  of  Figure  4.5  (b) . Figure  4.15(a) 
shows  the  resultant  loading/unloading  domains  and  Figure  4.15  (b)  the 
e-t  and  Ay  - t plots  for  x = 1.88  inches  and  x = 5.75  inches.  Com- 
paring Figure  4.15  (b)  with  Figure  4.9  (the  only  difference  between 

% 

the  two  problems  they  represent  being  in  the  presence  of  a wall  con- 
strained to  rotation  in  the  case  of  Figure  4.15  (b)),  it  is  seen  that 
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T (psi) 


(a)  Stress  trajectory  at  station  I (x  = 1.88") 


T (psi)  Final  yield  surface 


-CT  (psi) 

(b)  Stress  trajectory  at  station  II  (x  = 5.75") 


Figure  4.14  Stress  trajectories  at  stations  I and  II 
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the  effect  of  the  wall  is  strongly  felt  on  the  e-response  of  the 
second  interior  station  which  is  near  this  wall.  Other  changes  in 
the  results  for  the  semiinfinite  tube  and  the  finite  tube  are  of  a 
much  less  degree,  at  least  up  to  the  time  of  the  end  of  computations. 


CHAPTER  5 


CONCLUDING  REMARKS 

Numerical  solutions  based  on  rate-dependent  constitutive  models 
have  been  presented  to  compare  with  published  experimental  data  on 
plastic  torsional  stress  waves  and  combined  longitudinal  and  torsional 
waves  in  prestressed  thin-walled  tubes.  Some  typical  unloading  problems 
in  combined  stress  waves  have  also  been  considered.  The  effect  of 
radial  inertia  has  been  neglected  in  all  cases. 

Results  from  the  analysis  of  torsional  waves  in  copper  show 
that  a rate-dependent  theory,  in  addition  to  correctly  predicting  the 
incremental  strain  wave  velocity,  actually  gives  a better  description 
of  the  propagation  of  large  level  strains  than  a rate-independent 
theory  based  on  the  quasistatic  stress-strain  behavior.  Two  semilinear 
models  incorporating  linear  and  exponential  overstress,  and  a quasilinear 
model  have  been  examined,  and,  while  all  the  three  give  reasonable 
agreement  with  the  experiment,  the  quasilinear  constitutive  equation 
gives  slightly  better  overall  results. 

The  study  on  combined  stress  loading  waves  in  aluminum  shows 
that  a dynamic  material  response  that  can  be  adequately  represented 
by  a rate-independent  theory  can  also  be  represented  as  well  by  a 
rate-dependent  theory  with  suitable  rate-sensitivity  parameters.  It 
has  been  shown,  in  fact,  that  the  rate-dependent  solution,  while 
retaining  the  good  features  of  the  rate-independent  solution,  gives 
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a better  representation  of  the  velocity  of  propagation  of  various 
levels  of  incremental  shear  strain.  On  certain  other  aspects  of  the 
experimental  results,  like  propagation  speed  of  longitudinal  strains, 
both  the  theories  disagree,  though  sometimes  in  different  directions. 

The  work  on  unloading  considerations  in  combined  stress  waves 
shows  that,  as  a numerical  exercise,  it  is  a much  less  difficult  prob- 
lem than  the  corresponding  problem  in  rate-independent  theory.  Also, 
because  of  the  relaxation  process  involved  before  elastic  unloading 
can  occur,  the  initial  slope  and  the  instant  of  initiation  of  the 
loading/unloading  boundary  is  usually  different  from  the  rate- 
independent  situation.  Numerical  solutions  of  the  particular  problems 
considered  on  unloading  in  semiinfinite  and  finite  tubes  have  brought 
out  some  interesting  details  that  are  not  otherwise  easy  to  visualize. 

The  principal  import  of  the  present  work  is  numerical.  In 
terms  of  mathematical  methods,  the  present  study  can  be  seen  to  be 
a numerical  analysis  of  some  nonlinear  stress  wave  propagation  problems 
The  integro-diff erential  scheme,  the  higher  order  Courant,  Isaacson, 

Rees  method  (with  more  iterations  found  necessary  than  originally  pro- 
posed) , and  the  five-dimensional  Newton-Raphson  scheme  using  Gauss- 
Jordan  elimination  are  the  three  striking  features  in  the  present 
analysis.  The  convergence  difficulty  faced  with  conventional  itera- 
tion schemes,  when  representing  an  essentially  rate-independent  behav- 
ior by  means  of  a rate-dependent  theory,  can  be  greatly  surmounted  by 

% 

the  use  of  these  powerful  techniques. 

The  primary  need  for  research  particularly  in  combined  stress 
waves,  at  this  stage,  seems  to  be  reliable  data — more  in  number,  and 
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measured  with  different  experimental  techniques  for  various  prestress 
conditions  and  materials.  On  the  theoretical  side  there  is  a need  to 
develop  more  realistic  constitutive  models  that  can  predict  better 
the  roles  of  different  plastic  strain-rate  components  that  participate 
in  the  complete  set  of  equations.  On  the  computational  side  there 
remains  always  the  challenge  of  efficiently  solving  equations  in  two 
or  three  space  dimensions. 


% 


APPENDIX 


A. 1 Recursive  Scheme  for  Evaluating  the  Integrals 
in  the  Integro-Dif f erential  Approach  Given 
in  Section  2.7 

We  consider  the  integral:  I = J exp  ^g(s)^  ds 

o 

and  the  resultant  function:  cp  = A Sm.  (a  + 3D 

where  the  meaning  of  g(s) , A,  a,  3,  as  written  here,  is  clear  from 
Section  2.7.  To  prevent  the  integrand  in  the  above  integral  from 
exceeding  the  largest  number  that  the  computer  can  handle  (exp  (174) 
for  IBM  360/65)  for  any  t < t (where  tf  is  the  final  time) , the 
following  artifice  was  used. 

I.  Choose  some  t for  which  g(t^)  < 174,  with  a little  margin. 

For  t < t < t where  (g(t  ) - g(t  ))  < 174,  the  integral  may 
be  rewritten  as 


r ± r%L 

It  = exp  (g(t]L))[  exp  (-g(t  1))f  exp  (g(s))ds+  J exp  (g(s) 


- g(tx)) 


ds^J. 


Define  R^:  It  = exp  (g^))  Rt  > ^ 

.'.  cp  = Ag(t1)  + A Bn  (o-  exp  (-g^))  + 3Rt  > t ) • 

II.  Choose  t such  that  (g(t  ) - g(t  ))  < 174.  Then  for 
*2  < * * *3 
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t 

I = exp  (g(t  ))R  + exp  (g(t  ))  f exp  (g(s)  - g(t  ))  ds 

* 

r t 

= exp  (g(t2))|^exp  (~g(t2)  + gd^))!^  + f exp  (g(s) 

2 t2 

- g(t2))  dsJ 

= exp  (g(t2)Rt>t 

2 


.'.  cp  = Ag(t2)  + A 0n  (a  exp  (-g(t2))  + PRt>t  )• 

2 

To  prevent  a large  negative  exponent,  we  may  replace  exp  (-g(t2))  by 
exp  (-gCt^))  without  any  loss  of  accuracy.  This  recursive  procedure 
was  carried  up  to  13  stages  in  the  actual  solution. 
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A.  2 Flow  Chart  for  the  Integro-Diff erential 
Approach  Discussed  in  Section  2.7 


Typically,  I = integral  (see  Appendix  A.l)  evaluated  at  point  M. 
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A.  3 Flow  Chart  for  the  Modified  Courant , 
Isaacson,  Rees  Method  (Section  2.8) 
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A.  4 Chart  for  Combined  Stress  Waves 
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Schematic  Structure  of  the  Boundary 
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Gauss-Jordan  Elimination  Subroutine 


Details  of  Newton -Raphson  Scheme 
Used  in  Boundary  and  Interior 


The  actual  Newton-Raphson  routine  is  that  shown  between  the  two 


dashed  lines. 
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